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ABSTRACT 

Baker, James Richard. Ph.D. , Purdue diversity, May 1975. 
GEOMETRIC ANALYSIS AND RESTITUTION OF DIGITAL MULTISPECTRAL SCANNER 
DATA ARRAYS. Major Professor: Edward M. Mikhail. 

This thesis contains the results of an investigation performed 
in order to define causes of geometric defects within digital multi- 
spectral scanner (MSS) data arrays, to analyze the nature of the 
resulting geometric errors, and to investigate restitution methods 
to correct or reduce such geometric errors. 

The thesis includes a review of digital multi spectral scanning 
systems, including the recording and digital form for computer 
aided analysis. Causes of geometric errors within the data arrays 
are presented, and previous .investigations ; of geometric aspects of 
remote sensing systems are reviewed. 

The introduction of geometric transformation relationships 
for scanned data, from which collinearity equations for MSS may he 
derived, serves as the basis of parametric methods of analysis and 
restitution of MSS digital data arrays. The linearization of these 
collinearity equations is presented, including consideration of 
the functional assumptions made in order to model the stochastic 
changes in the exterior orientation of the sensor down the flight 


line . 
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Parametric algorithms for analysis end restitution based upon 
the above analytical treatment which were considered are: the 

direct use of MSS collinearity equations, and the use of piecewise 
polynomials based upon the linearized collinearity equations. 

In addition, nonparametric algorithms for restitution are introduced, 
as an alternative to the parametric algorithms. 

A proposed system for the geometric analysis and restitution 
of MSS digital data arrays is introduced. This procedure was used 
to test the methods of analysis and restitution, utilizing actual 
MSS dat# arrays from two aircraft flights . The results of these 
tests indicate that the collinearity equations can yield acceptable 
results when utilized for the analysis and restitution of such 
arrays. The investigation indicates that the piecewise polynomial 
algorithms are in general inferior to other methods of restitution. 

The nonparametric algorithms show great promise for the restitution 
of these arrays since the resulting accuracy of restitution is 
comparable to or better than that using collinearity, and they posses* 
definite advantages in computational efficiency. The arithmetic 
mean algorithm appears to perhaps represent a particularly efficient 
nonparametric algorithm for restitution of MSS digital data, at 
least for the arrays tested. 
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INTRODUCTION 

This introduction serves to define the problem to be investigated. 
The objectives of the research are presented, and the scope of the 
investigation is defined. An outline of the thesis is included as 
an aid to the reader. 

Statement of the Problem 

Computer aided analysis of digitally recorded multispectral 
scanner (I4SS) data arrays with the aid of high speed digital computers 
has made great progress recently. However, the digital data arrays 
which serve as the basis for these investigations are subject to 
geometric distortions, which are presently ignored or subordinated 
to the interpretation aspects. Consequently, the research problem 
of this thesis is to investigate the geometric distortions present 
within digital MBS arrays, and to formulate procedures for reducing 
the geometric errors present within the data. 

Research Objectives 

Based upon the problem statement above, the following research 
objectives were formulated. 

1.) To systematically review the methods of recording MSS data, 
generating digital MSS data arrays, and performing computer 
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aided analysis* in order to define the context within which 
the geometric aspects of these data are considered. 

2. ) To define the geometric distortions present within the 

data arrays. 

3. ) To analyze the transformation equations of MSS data recording 

and derive from them functional forms for investigating 
geometric errors within these data and algorithms for re- 
ducing these errors . 

4. ) To introduce nonparametric methods not based upon these 

transformation relationships as an alternate method of 
reducing geometric error within the arrays . 

5. ) To identify the problems associated with these methods 

when working with digital rather than continuously recorded 
information and to present methods to deal with these problems. 

6. ) To test the analysis methods of 3. and 4. above on real 

MSS digital data arrays. 

Scope of the Investigation 

To satisfy the research objectives above, the investigation 
is limited to evaluation of MSS digital data arrays obtained by 
airborne multispectral scanners. However many of the concepts and 
procedures presented will be applicable to data gathered from space- 
craft. The investigation was limited to actual data arrays presently 
existing. No simulation of data was attempted. The investigation 
includes both theoretical modeling and numerical analysis by computer 
as an integral part of the research. 
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Outline of the Thesis 

The thesis is organized in such a manner as to provide a logical 
sequence, in general progressing from the more theoretical aspects 
of analysis to the more practical applications utilizing the relation- 
ships introduced. In Chapter 1 a "brief introduction to the field 
of remote sensing, the most widely used remote sensing systems, and 
basic radiometric concepts common to many of them is given. The 
physical configuration of the multispectral scanner and the processing 
of the information gathered by the scanner into digital data arrays 
are presented. Techniques of computer aided analysis utilizing the 
digital arrays thus generated are introduced. Also included in this 
chapter are derivations of the basic sources and types of geometric 
errors which have been included in the arrays. The chapter concludes 
with a review of previous investigations of the geometric aspects 
of remote sensing systems. 

Chapter 2 begins with the introduction of a general remote 
sensing transformation, including both radiometric and geometric 
information. The geometric portion of this transformation is derived 
in its most general form for the MSS system. From this basic geo- 
metric transformation, the collinearity equations for multispectral 
scanners are derived. The procedure for linearizing these equations 
by Taylor's series for performing resection by least squares is 
presented. This linearization includes the dynamic effects of time 
varying orientation parameters. 

In Chapter 3 restitution methods for MSS data arrays are discussed. 
The purpose of these restitution methods is to process the data to 


form arrays which are essentially orthographic projections. First, 
simple resampling algorithms are presented which are useful as 
pre-processing steps. The direct use of MSS collinearity equations 
as a restitution method is next presented. The lose of polynomial 
forms for restitution is given, and a discussion follows of some of 
the problems associated with using piecewise polynomials and the 
collinearity equations. A completely different approach to restitution 
of MBS digital data arrays is then introduced, that of nonparametric 
methods. The basic forms for these nonparametric restitution al- 
gorithms are presented along with the major advantages and disad- 
vantages of each. The chapter next deals with the unique problems 
associated with transforming digital data arrays into the desired 
form. The chapter concludes with a discussion of the problems 
associated with determining elevations within digital data arrays 
which have been obtained from singly scanned areas. 

Chapter 4 deals with the testing of the above methods of analysis 
and restitution on actual MSS data arrays . Two flight lines of data 
were obtained and analyzed using each of the above methods. The 
results are statistically analyzed and compared. 

The thesis concludes in Chapter 5 by summarizing the research 
effort, stating the conclusions drawn from the investigations performed, 
and enumerating recommendations for further research in this area. 
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1. THE MULTISPECTRAL SCANNING (MSS) SYSTEM 

This chapter is intended to provide the background material 
required for an understanding of the investigation. A brief intro- 
duction to the discipline of remote sensing, and a brief survey of 
remote sensing systems is provided. The basic radiometric concepts 
required for an understanding of recording data by multispectral 
scanners is described, and the physical configuration of the scanner 
is presented. The digital aspects of multispectral scanning (MSS) 
systems are presented, including a brief introduction to the field 
of machine aided interpretation utilizing digitally recorded data. 
The basic geometric distortions inherent in the digital arrays are 
presented in order to obtain an appreciation of geometric defects 
within the data. The chapter concludes with a review of previous 
investigations in geometric analysis of remote sensing systems. 

1.1 Introductory Remarks 

The term "remote sensing" may be applied to virtually any 
activity in which information is gathered without physical contact. 
As generally used, however, the term is usually applied to those 
systems which yield data about the Earth’s surface and its resources 
from a considerable distance. References [l] through [4] contain 
definitions of remote sensing systems by persons prominent in this 
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field. The definition from [^] is one of the most descriptive: 

"Remote sensing denotes the joint effects of employing modern sensors, 
data processing equipment, information theory and processing method- 
ology* communications theory and devices, space and airborne vehicles, 
and large systems theory and practice for the purpose of carrying out 
aerial or space surveys of the Earth's surface." It may be seen 
from this definition that the existing disciplines of photogrammetry, 
radar gramme try, hologrammetry , and many other activities can be 
considered as divisions within remote sensing. 

Remote sensing systems may be classified according to many 
different criteria [ 5 ]: geometric properties, stationary or moving 

platform, image forming or non image forming, and passive or active 
with respect to the source of a carrier field such as electromagnetic 
energy, sound, or gravity. Figure 1.1 depicts the various classifi- 
cation schemes, ending with examples of sensors which fall into the 
various categories. 

An active-passive differentiation between sensor types depends 
upon whether the sensor generates its own energy source for the 
carrier field (active) or relies upon an external source for the 
carrier field (passive). Both active and passive systems may have 
image forming capabilities. Classification as regards geometry 
would include area recording (e.g. frame camera), line recording 
(e.g. panoramic and continuous strip cameras), or point recording 
(e.g. optical mechanical scanners). References [6], [T]» and [8] 
contain other classification schemes or tabular data concerning the 
properties of remote sensors. 


Remote Sensing Systems 
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Figure 1.1. Remote Sensing Systems 
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The subject of this investigation, the multispectral scanner, 
is a passive system operating from a moving platform (either aircraft 
or spacecraft) with potential image formation and records in a point 
by point mode. The term "point" is used here to mean a small but 
finite area representing the system's resolution element as will 
be explained in detail later in this chapter. Scanning modes for 
ICS systems may be linear about a near vertical axis, linear about 
an axis tilted with respect to vertical, or conical about an axis 
which may be vertical or tilted. This investigation will be limited 
to the case of linear scanning about a near vertical axis in a 
direction normal to that of the motion of the platform. 

The carrier field for the ICS system is composed of electro- 
magnetic waves. Figure 1.2 presents the spectrum of the electro- 
magnetic fields . The human eye is sensitive to the visible portion 
of the spectrum, between the wavelengths of O.U to 0.7 micrometers. 
The entire band between 0.3 and 15 micrometers is referred to as 
the optical portion of the electromagnetic spectrum, since electro- 
magnetic waves may be refracted and focused using lenses and prisms 
in this region. It is within this range that most passive remote 
sensing systems operate. The region of the spectrum between 0.3 
and 1.0 micrometers is called the photographic region where photo- 
graphic emulsions may be used to record data. Wavelengths below 
0.^ micrometers are in the ultraviolet region of the spectrum, and 
that portion above 0.7 micrometers is the infrared portion. The 
region between 0.7 and 3 micrometers is referred to as the reflective 
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infrared region, and the region from 3 to 15 micrometers is called 
the emissive region, since energy in these wavelengths is dominated 
by emission from a body as a result of thermal activity at temperatures 
of approximately 300° K, rather than being reflected. 

For the subsequent discussion of basic radiometric relationships, 
the term power is used to denote the time rate of energy emanating 
from an object, or incident upon some sensor. Radiation power would 
then denote the time rate of energy due to reflected and/or emitted 
electromagnetic energy. Radiance may be defined as the radiation 
power per unit area directed into a unit solid angle from an object. 
Irradiance is used to denote the radiant power per unit area incident 
upon an object or sensor. Reflectance is a unitless quantity between 
zero and one and represents the portion of energy incident upon an 
object which is reflected from it. 

The source of the energy incident upon the sensor's detectors 
is reflected and emitted radiant energy from the scene being sensed. 

All energy coming to Earth from the sun is either reflected, scattered 
or absorbed and subsequently emitted by objects on Earth. The total 
radiance from an object is composed of two components, reflected 
radiance and emitted radiance. In general, it may be stated that 
the reflected radiance forms a dominant portion of the total radiance 
from an object at shorter wavelengths of the electromagnetic spectrum, 
while the emissive radiance becomes greater at longer wavelengths [9» 
Figure 1.3 depicts in schematic form the interaction between the 
radiant energy emanating from the sun and a ground scene. In this 
figure, E denotes emitted radiation, and R denotes reflected radiation. 
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The spectral radiance, in watts/sq. cm-steradian-micrometer 
from a scene in the reflectance region may be given as [9] 

= ]L E^ p cos 0^ (1.1) 

TT 

in which E^ is the solar irradiance in watts/sq. cm incident on the 
object at a wavelength X, 0 g the sun angle and p the reflectance 
of the scene, a unitless quantity. Figure l.k shows a plot of ra- 
diance versus wavelength for typical agricultural scene, in the 
reflectance region, having wavelengths between 0.3 and 3.0 micro- 
meters [9]. The figure shows the reflected radiance as it would 
be recorded near ground level. For this plot, p is assumed equal 
to 0.1, a reasonable assumption for vegetation in the visible portion 
of the electromagnetic spectrum or for soils in the near infrared 
portion of the spectrum. 

In the emissive region of the spectrum, the radiance from an 
object may be represented by Planck’s Lew, 

L, = (1.19 x 10 4 )e (1.2) 

X 5 [exp ( 14338/XT )-l ) 

in which T represents the absolute temperature of the object in 
degrees Kelvin, and e represents the emissivity, a unitless coefficient 
which is unity at all wavelengths for objects called black bodies. 
Objects are called gray bodies if their emissivity is a constant 
between zero and one. Most objects however, are neither black bodies 
nor gray bodies, and the emittance may be a complicated function 
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of wavelength. Figure 1.5 shows a plot of spectral radiance for 
a typical black body, emittarice equal 1, at an absolute temperature 
of 300° K. 

The combination of two curves from the reflectance and emittance 
regions, then, would represent the total spectral response of the 
object. Figure 1.6 shows a simplified example in which the agri- 
cultural scene represented in the reflective region by Figure 1.4, 
is assumed to behave as a black body in the emissive region, represented 
in Figure 1.5. It is this plot which represents the "spectral 
signature" of an object, and it is the di^.'erence between such 
signatures for different objects which allows the differentiation 
between objects using the multispectral approach of remote sensing. 

The summation of the preceding plots will yield the magnitude 
of radiance emanating from an object. In order to evaluate the 
energy* or power (time rate of energy) incident at the sensor, how- 
ever, the transmission through the atmosphere must also be considered. 
The atmosphere will attenuate and modify energy very significantly 
by molecular absorption for many spectral wavelengths. There are, 
however, spectral intervals at which the atmoshpere is reasonably 
"transparent" to electromagnetic radiation. These transparent spectral 
regions are called "windows", and sensing from a remote location 
requires that these spectral intervals be utilized. Figure 1.7 
presents a graph of the transmission of the atmosphere for the 
wavelengths of the optical portion of the spectrum. The plot does 
not include the effects of atmospheric scattering. Table 1.1 sum- 
marizes some of the more important atmospheric windows which have 



Spectral 
L x (watts /cm 2 - nr 



Figure 1.5. Spectral Radiance in the Emissive Region [ 9 ] 




Wavelength (micrometers) 

Figure 1.6. Total Spectral Radiance [ 9 ] 


Spectral Radiance, 

L ^ ( watts /cm - micrometer-steradian ) 
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Figure 1 . 7 . Atmospheric Transmission [11 3 
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been utilized in MSS remote sensing systems. The most often used 
windows in remote sensing are numbers 1, 2, 3, 4, 5, 6, and 9 [ill - 

Table 1.1 

Atmospheric Windows [ll] 


Window No. 

Wavelength Limits 

1 

.40 to 

.72 pm 

2 

.72 to 

.94 pm 

3 

.94 to 

1.13 pm 

4 

1.13 to 

1.38 pm 

5 

1.38 to 

1.90 pm 

6 

1.90 to 

2.70 pm 

7 

2.70 to 

4.30 pm 

8 

4.30 to 

6.00 pm 

9 

6.00 to 

15.00 pm 

10 

15.00 to 

25.00 pm 


If the atmospheric transmission is known for a spectral wave- 
length interval A to (A + A A) , the resultant power flow into the 
sensor may be given by [ 9 ] 

P =* x a L x A g 2 (AA) (1.3) 

in which P is the power in watts, t_ the transmission of the atmos- 
phere, and 3 the resolution angle of the system, defined as the 
angle from the sensor within which the electromagnetic energy is 
received at any instant, and A is the effective area of the receiving 


optics . 
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With these preliminary remarks in mind, attention is now directed 
to the multispectral scanner itself. In the following section a 
discussion is given of the manner in which these radiometric proper- 
ties of matter may be utilized to record analogue data using the 
multispectral scanner, and how these analogue records may be processed 
to produce digital data arrays for analysis. 

1.2 Digital Multispectral Scanning System 
In this section, the steps taken to record the electromagnetic 
energy using a multispectral scanner will be discussed. The processing 
steps which take the data from electromagnetic waves , to analogue 
records from the MSS, and from these analogue records to digital 
data arrays will be discussed. The final form of data, the digital 
array is the form which is utilized at the Laboratory for Applications 
of Remote Sensing (LARS) to perform automated interpretation, and 
it was this data format which was used to investigate the geometric 
properties of the digital MSS system. 

1.2.1 The Multispectral Scanner 
Figure 1.8 depicts the scanning and data collection system 
for a typical multi-spectral scanning system. The system depicted 
is an aircraft system, in which a rotating mirror or prism scans 
the terrain in narrow strips or lines oriented normal to the direction 
of flight. Forward motion of the aircraft provides continuous 
coverage by assuring an advance for each line, allowing for some 
overlap between successive strips. The system utilized at LARS 
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to generate basic data for digital processing is discussed 
further. 

The electromagnetic energy incident upon the mirror or prism 
is optically focused upon a separating prism or diffraction grating 
which separates the incident energy into different wavelength bands. 

Each band is directed to an appropriate detector, and the signal 
generated is recorded on magnetic tape as an analogue signal, each 
spectral band representing a data "channel". The data is recorded 
on l4 track one inch magnetic tape as a series of continuous signals. 

The first track of tape contains a synchronizing signal to note the 
beginning and end of each scanner revolution. Also on this track 
is a signal from a roll stabilization gyro, which indicates the true 
roll of the aircraft during each scan line period. The remaining 
13 tracks are data channels, which record calibration sources both 
internal and external to the aircraft , as well as recording the 
radiance from the ground scene of interest. Therefore, up to 
13 channels in 13 spectral bands may be recorded simultaneously. 

The optical system shown in Figure 1.8 is of the converging optics 
type coupled with a rotating mirror-prism. Other optical configura- 
tions are available utilizing an oblique mirror or wedge, and 
converging optics [10, 11] • If digital analysis of MSS data is desired, 
these tapes are subject to digitization in an analogue to digital 


sampling procedure. 
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1.2.2 Advantages of Multispectral Recording 

The M3S sensing system has several advantages over more con- 
ventional sensors such as photography. The data are recorded in 
a quantitative manner suitable for processing into numerical data. 

The recording and processing can be calibrated, resulting in high « 

repeatability of decisions based upon statistical procedures. 

The possibility of detailed numerical analysis allows the discrimination 
of subtle differences within the data resulting in a high probability 
of isolating classes of interest through automated interpretation 
and classification. A class may be thought of as any physical divi- 
sion of interest to the investigator, which is believed to be spectrally 
separable from other classes. Examples of classes are wheat, com, 
oats, and other crops which may be appropriate in an agricultural 
experiment utilizing airborne sensors. Another example of classes 
would be bare soil, green vegetation, and water, which may be appro- 
priate for an investigation of a rural environment using spacecraft 
gathered data. Urban, suburban, and rural may define classes of 
interest for studies of urban areas using spacecraft data. 

An illustrative example is cited by Holter [4], in which a 
comparison of the information available on black and white film, 
color film, and the multispectral system is made. For this example 
it is assumed that the electromagnetic spectrum is divided into 
12 wavelength bands, and each band is assumed to have 10 levels of 
intensity, or density. Black and white photography senses 
instantaneously all of the bands in the photographic region of 
the electromagnetic spectrum. The maximum density would occur when 
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the maximum intensity is present in each band, and there will be 

10 x 12, or 120 density values discernible below this value. Each 

of these values , however may occur from a different combination 

of bands , and the band combination causing this density cannot be 

discerned. Color film may be thought of as recording in 3 bands, 

representing the 3 emulsions used. If each emulsion is assumed 

sensitive to one third of the 12 spectral bands , then each emulsion 

could discern 1+0 densities, and the discernible number of states 

rises to approximately 6.4 x 10 . A multi spectral sensor,' on the 

other hand, senses each band separately, and stores the data separately, 

so that the total number of states available from the multispectral 

12 

sensor data which could be recognized is 10 

Although meaningful interpretation can and has been done by 
conventional interpretation methods utilizing multispectral images 
in order to make efficient use of all of the data available in this 
form, human decision making becomes too slow by several orders of 
magnitude. Further, the human interpreter finds it impossible to 
simultaneously discern differences in spectral tones, as displayed 
upon images, for more than a few (3) spectral bands. Such data may 
be analyzed with automated systems of interpretation, utilizing 
digital or analogue computing equipment. If the digital approach 
is chosen, the analogue flight record must be converted into data 
forms compatible with the digital computer, such as digital data 
tapes. If such a form is generated, and appropriate digital decision 
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algorithms are written, then computer aided analysis may he 
performed with speed and accuracy, and quantitative results 

obtained . v 

1.2.3 Generating Digital Data Arrays 
In order to utilize digital data handling and interpretation 
of MSS data, a system for converting the analogue record to digital j 

data format must he imp lamented. One such system, flown by the 
Environmental Research Institute of Michigan (ERIM) relies upon 
a roll gyro signal to establish the roll of the aircraft for each 
scan line. Figure 1.9 (adapted from [12]) depicts the manner in 
which the roll of the aircraft is monitored in flight and later 
utilized in analogue-digital (A-D) conversion to generate roll com- 
pensated data. As the scanning prism rotates, a magnet attached !„ 

to the shaft passes a coil at the same point for each revolution. 

Referring to Figure 1.9a, as the magnet passes the coil, it generates 
a synchronizing pulse which is placed on the analogue record (Fig- 
ure 1.9b) , and simultaneously triggers a constant slope electronic 
ramp signal generator. When the signal generated is equal to that 
of the output signal from the roll gyro (Figure 1.9c), indicating 
the position of the aircraft roll with respect to nadir, a roll pulse 
is generated with respect to the synchronizing pulse upon the analogue 
record. Thus the magnitude of roll is recorded upon the analogue 
record for each scan line. These synchronizing and roll pulses are 
recorded upon the first track of the analogue recording tape within I 


the aircraft. 





Figure 1.9. Digital Sampling of Analogue Data [12] 
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As the scan continues past the aircraft aperture, the spectral 
radiance is recorded across the ground scene on the remaining tracks 
of the tape for up to 13 data channels. As the rotation is continued, 
the sensor receives reference signals (from a light source, heat 
plates for the emissive regions of the infrared spectrum, and the 
sun). These reference values are used subsequently as calibration 
data. Figure 1.9b depicts the first track and two typical data 
channels . 

Digital sampling is then done in the laboratory using analogue 
to digital (A-D) conversion equipment. The size of a data sample, 
and hence the angular interval of a digitized element in the along- 
scan direction is determined from the angular scan rate of the rotating 
prism, by selecting an appropriate sampling rate in the A-D converter. 
As an example, consider a data set in which the scene was scanned 
at a rate of 3000 revolutions per minute. The resulting period 
of rotation for the scanning prism, t r , is then 1/50 sec. , and 
is represented on the analogue record by the distance between sucessive 
synchronizing pulses, as shown in Figure 1.9b. If a nominal desired 
angular sampling interval y, is desired, the number of samples between 
synchronizing pulses may be computed and a time interval of sampling 
computed, which may be dialed into the A-D converter. 

For a nominal angular sample interval, or instantaneous field 
of view, of y - 6 milliradians, the number of samples between syn- 
chronizing pulses (representing one revolution, or 2 ir radians) 
may be computed as 

n - 2 tt = 1048 samples, 

.006 
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or, in terms of time interval of sampling in the A-D converter. 

At = .006 (.02) = 1.9 x 1CT 5 sec. 

2 n 

In practice, the time interval of sampling is limited to stepped 
values, and hence the actual sampling resolution will he slightly 
different from the nominal value. 

After the time interval for sampling is known, the actual 
sampling is done with respect to the roll pulse. Beginning at the 
roll pulse position, the analogue record for each channel is instan- 
taneously sampled at the time intervals calculated, as shown in 
Figure 1.9t>* Since the sampling begins at the roll pulse position 
for each scan line, any sample, say sample number 100, on each scan 
line is recorded on the digital tape such that it represents the 
same angle with respect to the nadir, resulting in roll stabilized 
data. Hie analogue data for each channel at every sample point 
is quantified on a scale ranging from 0 to 255 > in order that it 
may be represented on the digital data tape (or in the digital computer) 
by a single byte ( 8 bits ) . A detailed explanation of the method 
in which the resulting data are stored on computer compatible tapes 
is given in reference [12], 

For each resolution element, there is a position in a matrix 
representing the scan line number in one dimension, and the element 
number within that line in the other dimension. Associated with 
that position are several spectral values equal in number to the 
spectral bands, or channels, of the system. The data, then, may 
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be considered as a three dimensional matrix of spectral values of 
size i x j x n, where i is the number of scan lines in the flight, 

J is the number of elements within each scan line , and n is the number 
of data channels recorded. Subsequent automated classification 
and interpretation schemes then consider an n-dimensional vector 
for each element in the i x j positional array, as the data form. 

In the following section, the use of the data arrays thus gen- 
erated to form images is discussed. The use of computer aided analysis 
algorithms operating on the data within these arrays is also presented. 
The development and application of these computer aided analysis 
and classification procedures constitutes the ma^or portion of the 
current research and development tasks performed at LARS. 

1.3 Computer Aided Analysis Using Digital MSS Arrays 

In the preceding section, the generation of digital MSS data 
arrays was discussed. In this section, the use of these arrays 
to extract meaningful information concerning resources at or near 
the Earth’s surface will be presented. 

1.3.1 Displays of MSS Data Arrays 

Whether the arrays generated are to be utilized for interpretive 
purposes or to investigate geometry, it is convenient to be able 
to display the information visually. The task of machine aided 
interpretation utilizes the arrays directly, after preliminary 
decisions are made from displays of the data. An alternate method 
of considering the data format, which is of more interest from the 
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photogrammetric aspect, is to consider the data as an image. Under 
ideal conditions, the image formed by a display of the digitized 
spectral values may be seen to consist of two types of projections 
as will be discussed in a later section of this chapter. In the 
direction of flight, the data is recorded in a line-by-line fashion, 
resulting in an essentially orthographic projection. In the direction 
perpendicular to the direction of flight, the image is a perspective 
projection about the effective perspective center of the scanner. 

If a display of this image is desired, the information stored 
on magnetic tape, must be processed to form a gray scale matrix. 

In the LARS system, the data analyst determines the number of gray 
scale levels desired (usually 10 to 12), and histograms are calculated 
for all spectral values within the area of interest. 

These histograms are then divided into a number of bins equal 
to the number of gray scales desired. The abscissa is divided such 
that the same number of data array elements fall into each bin. 

Each interval along the x axis is then assigned a symbol or gray 
scale level. When the array is subsequently displayed, the result 
is that every symbol will occur with the same frequency within the 
display. The result is equally active gray scale levels over the 
image. Display of the image may be accomplished in two ways. 

The data tapes may be used to drive the scanning raster upon a 
television monitor, with the gray scale levels determining 
the voltage of the raster signal. This form of display utilizes the 
LARS Digital Video Display Unit, shown in Figure 1.10 [13], and 
henceforth referred to as the digital display unit. An alternate 
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display form is a gray scale printout utilizing the line printer 
of the computer. Different gray scale values are represented by 
alphanumeric characters available on the line printer. Figure 1.11 [ 13 ] 
illustrates these two methods of display. The digital display unit 
is well suited to the display of gray tone imagery and allows great 
flexibility, through the use of an attached light pen, for purposes 
of interpretation. The line printer image form is useful for geometric 
studies, since it readily allows for the isolation of single element 
positions. 


1.3.2 Computer Aided Analysis 

As mentioned previously, the abundance of data generated by 
MSS systems require that interpretation be automated as fully as 
possible , to make maximum use of the unique properties of the multi- 
spectral approach. Computer aided analysis in this context, 
will denote the assignment of each element in the data array into 
one of a finite number of pre-assigned "classes" of interest, according 
to some decision rule. Different algorithms for this purpose abound 
in literature [lU through 193* A brief introduction to the method 
used at LARS will be presented in this section. The method is based 
upon statistical pattern recognition techniques. 

Figure 1.12 [ lU ] depicts in a block diagram a general pattern 
recognition system. The receptor or sensor in this case would be 
the multispectral scanning apparatus . Data from this sensor may be 
represented as an n-dimensional vector for each element in the posi- 
tional array, where n would be the number of spectral channels used. 
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Figure 1.12. A Pattern Recognition System [l!+] 
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The purpose of the feature extractor is to generate a feature space 
vector containing maximal information for discrimination purposes. 

As an example which has been used in the LARS system, the data analyst 
selects the desired dimension of the feature space vectors to be 
used. The analyst may elect to perform classifications using four 
of the channels when originally a total of twelve channels were 
recorded. A statistical quantity called "divergence” between class 
pairs is computed for all possible groups of four channels included 
within the original twelve channels. This divergence is a measure 
of the spectral separability of the classes based upon each group 
of four channels. Based upon this quantity, one group of four 
channels would be selected to perform the classification. Usually 
the four channel combination chosen will be the set which has the 
greatest average divergence between pairs, although other selection 
criteria may be used. The net result of the feature extractor as 
used at LARS is to select some subset of dimension m from the original 
n channels of the measurement space vector (m ^ n). The feature 
extractor could, however, be a more complicated transformation. 

The measurement space vector would be 
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and the feature space vector for analysis would be 



(1.5) 


The decision maker is of course an important step in the 
classifying procedure, since it assigns each point in the data array 
to some particular class based upon a decision rule formulated from 
calculations made upon the feature space vector. 

Each element and its associated m-dimensional vector represents 
a point in an m-dimensional feature space, and it is the location 
of points within this space which serves as a basis for classification. 
Figures 1.13 and l.lU represent a two dimensional example, in which 
two wave length bands, X ls and are the feature space coordinates. 
Figure 1.13 (course notes, EE 595, Purdue University) shows typical 
spectra of soil and vegetation. The spectra are assumed to be 
sampled as noted on the figure at wavelengths X-j_ and X 2 . The points 
in the two dimensional feature space (X^_, Xg) are then plotted in 
a two dimensional feature space as shown in Figure l.lU. 

Natural phenomena, however, such as radiated and reflected 
electromagnetic energy from which spectral data are derived exhibit 
some inherent randomness. Spectra from two separate representatives 
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Figure 1.13. Spectral Reflectance for Vegetation and Soil 

(From course notes, E.E. 595 » Purdue University) 
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Figure l.lU. Representation of Two Classes in a 
Two Dimensional Feature Space [iL] 
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of the same class would not be identical, and therefore the X^) 

measurements would not fall at the same point in the feature space. 
Instead, samples from different elements tend to fall in areas 
or "clouds" within the feature space. For example, Figure 1.15 
depicts many observations from 3 classes plotted in the feature 
space. The decision making function of the classifier divides the 
feature space into a number of regions (the basis of which will be 
discussed subsequently) each representing a class of interest. 

The decision making algorithm examines each new data point to determine 
into which region it will fall, and classifies the point accordingly. 

The classification problem then, becomes one of defining decision 
surfaces separating the appropriate regions, as shown in Figure 1.15. 
For the two dimensional feature space shown, the decision surfaces 
become lines. In order that the classifier may designate the decision 
surfaces before classification of unknown data points takes place, 
definition of these surfaces is done based upon "training samples", 
which may be generated from two different approaches, the data bank 
approach and the extrapolation approach [l]. 

In the data bank approach, an attempt is made to gather as 
many samples of spectral signatures as possible , in as many classes 
as may be useful to the data analyst. For any given problem a subset 
of classes present in the scene scanned would be chosen and decision 
surfaces calculated between these classes. Each incoming data 
point would then be placed into one of the regions defined by these 
classes. This method would have the advantage of requiring a minimum 
of a, priori knowledge of the data to be classified, and the same 
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Figure 1.15. Classification Regions for Three Classes 
in a Two Dimensional Feature Space [lU ] 
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data bank of spectral signatures could be used to classify any number 
of nev data sets . 

This method has serious drawbacks, however, since the data 
bank of spectral signatures would have to be very large, and would 
have to include signatures gathered under many different conditions 
in order to account for natural statistical variations, such as 
temporal (time) variations, sun angle variations, variations in 
the atmosphere, seasonal variations, etc. 

Using the extrapolation approach, however, these problems 
are circumvented to some extent . The training samples for this 
method are samples of known classification, drawn as a subset of 
the data to be classified, assuring that the training samples were 
gathered under the same conditions as the data to be classified. 
Decision boundaries are determined from these training samples, and 
the entire set of data is then classified accordingly. Thus, the 
classification is based upon an extrapolation of a few known points 
within the data set itself. 

The method has the advantage that much less calibration is 
required of the sensor. However, it does require the acquisition 
of some a priori knowledge of the scene being scanned. This knowledge 
may be gained from ground observations in populated or accessable 
areas, limited photographic missions in inaccessable areas, or a 
combination of the two. It is this method which is primarily utilized 
at LARS . 

The actual assignment of decision surfaces is the subject of 
continuing investigations. As a simple example, the centroid of each 
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class may be calculated, and the surface equidistant between centroids 
of adjacent classes used as the decision surface. In this case, 
the decision surfaces shown in Figure 1.15 would become straight 
line segments. At LARS, the decision surfaces are calculated to 
minimize the average expected loss of classification accuracy, as 
based upon some assumed loss function. The data in the spectral 
channels is assumed Gaussian in its randomness, and the resulting 
classification system is a maximum likelihood decision rule. 

If the feature space is divided into regions as Just described, 
every incoming point must be assigned to one of the designated 
classes. In virtually any scene, however, there are elements which 
do not belong reasonably to any of the finite number of classes 
designated. To avoid these gross mis classification errors, the 
concept of thresholding is used. If the probability of an element 
belonging to any of the designated classes is smaller than some 
preassigned threshold probability, then that element is assigned 
to none of the classes of interest but instead is placed in a null class 
which will represent "all others". 

If training sample areas are difficult to delineate, another 
procedure called clustering, or unsupervised classification, may 
be used to supplement the above procedure. In using clustering, 
the analyst decides upon some number of groups into which he desires 
that the spectral data be divided. The clustering algorithm then 
performs this division, dividing the data into this number of groups, 
choosing those groups which are most separable based upon the spectral 
values stored in the data arrays. These groups are called "clusters". 
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After these groups have been isolated, however, the data analyst 
must determine whether each cluster represents a class in the physical 
sense, whether each cluster represents a single class or a combination 
of classes, and whether the classes isolated are of interest. 

The above brief summary of data classification has used as 
an example the system utilized at LARS. The system is based upon 
spectral separability of data elements. It si Jcu-l be noted that 
other information may also be used for interpretation, such as spatial 
information, temporal (time) information, polarization, and combinations 
of these. 

To date, all of the algorithms and procedures described have 
used as the basic data source the arrays directly as they have come 
from the A-D conversion. No attempt has been made to analyze the 
geometry of the data arrays . This analysis is desirable in order 
to be able to extract from digital MSS data arrays not only the 
interpretive information of ’’what" is included in the information, 
but also the metric information of "where" and the quantitative 
information of "how much". In the next section an introduction to 
the geometric distortions present in the data arrays will be given, 
to aid in -understanding the types of geometric deformations to which 
the data arrays are subject. 

1.4 Geometric Distortions in MSS Imagery 

Investigations concerning ICS digital data arrays have previously 
dealt primarily with the automated interpretation phases of analysis 
through pattern recognition using statistical concepts. Geometry 
of the image has been either neglected entirely or subordinated to 
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the interpretation analysis. The MSS data, however, is subject 
to geometric distortions, which may be serious if it is desired 
to analyze the data to extract metric as well as qualitative information. 
It would be desirable to be able to obtain metric information from 
the data arrays. Positional information may then be associated 
with each properly classified element, and reliable information 
concerning areas may be obtained. 

1.4.1 Ideal Geometry 

Figures 1.16 and 1.17 illustrate geometry of a multispectral 
scanner under ideal conditions. The case illustrated represents 
the recording of the j-th element in an arbitrary scan line, i. 

In this case the following assumptions are made: 

1. ) The aircraft is flying perfectly straight, at a constant 

elevation above datum, and at a constant ground speed. 

2. ) The aircraft is subject to no angular exterior orientation 

perturbations. That is, there is assumed to be no roll (u>), 
pitch (<j>), or yaw (k) of the aircraft. 

3. ) Each scan line is assumed to be instantaneously recorded. 

Under these conditions the resulting image is recorded as an ortho- 
graphic projection in the direction of flight (x), and as a perspective 
projection in the direction normal to the direction of flight (Y). 

In Figure 1.16, 3 represents the angular resolution of the 
beam in the down strip or X direction and will be the same as the 
physical resolution of the scanner. The angle y represents the 
effective angular resolution in the along-scan or Y direction and 
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is determined during the analogue to digital conversion. This represents 
the effective instantaneous field of view (IFOV) of the sensor. 

The total angle scanned, called the total field of view (TFOV) is 2a. 

The effective aperture of the sensor is represented by the 3 by y 
solid angle. 

Figure 1.17a depicts the down strip or X direction. For this 
direction the Xj coordinate for the object point j is a function 
of time down the strip and may be written 

T i 

X- = x c = X Q + / V(t)dt (1.6) 

T o 

in which X c is the X coordinate of the scanner at the instant of 
imaging point j on scan line i, X^ for point j, T^_ is the time 
epoch at the instant of recording point j and is the same for all 
points on scan line i under the assumption that the scan line is 
instantaneously recorded, X Q is the X position of the first scan 
recorded and V is the aircraft velocity which is assumed constant 
for the ideal case. If the scanning rate is adjusted such that an 
overlap (overscan) will occur between scans at the nominal aircraft 
velocity, then the average advance for each scan line will be (l-S)dX, 
where S is some overlap factor, and dX is the scan width on the 
datum given by (see Figure 1.17) 

* Z c 3 (1.7) 


dX 
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In this expression Z Q represents the altitude of the aircraft above 
datum when recording scan line i. Further, the image element size 
in the down strip direction is given by 

dx = c3 (1.8) 

in which c is an effective sensor constant, analogous to the record- 
ing barrel radius (camera focal length) in panoramic photography. 

Noting from Figure 1.17a that 

V dt - (l-S)dX (1.9) 

and solving Equation (1.8) for 8, the following expression results 
from substitution of the expression for 3 into Equation (1.7). 

x 

Xj = X c = X Q + / (1-S) Zjvdx' (1.10) 

o c 

For T Q = 0, the image strip coordinate is equal to zero, and at 
time t = Tj_, x’ = x, the image x coordinate of point j. 

Applying the ideal assumptions mentioned above, in which aircraft 
velocity and flying height above datum are assumed constant down 
the flight line, the expression becomes 

Xj = X c = X c + (l-S) Zg_ x (1.11) 

c 

In reality, of course, the aircraft velocity and flying height 
may be a function of time, and hence of the image x coordinates. 

If a polynomial form for these quantities is assumed, then a polynomial 
form of one degree higher will result for X,. after integration of 
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Equation (1.10). If a harmonic form is assumed, then a harmonic 
form will remain after integration. 

In Figure 1.17b is shown the recording geometry in the direction 
perpendicular to the direction of flight. For this case, the ground 
coordinate Yj is given by 


Y 0 + *'3 
or 

Yj = Y c + (Z c - Z) tan6j 


( 1 . 12 ) 


in which Y^ is the ground coordinate of point j, Y c is the Y 
coordinate of the sensor at the instant of recording j , Z is the 
elevation of the terrain (assumed level) in which Zj in Figure 1.17b 
has been replaced by the constant elevation Z. The angle 0^ is the 
scan angle at the instant of recording given by 


0j * yj/c (1.13) 

and yj is the image position normal to the flight line axis and 
would be represented in digital recording by the column position 
with respect to the scan line center. 

1.4.2 Variations from Ideal Geometry 
This idealized geometry is obviously not realized in practice, 
since the ideal positions as given by the projection Equations (l.ll, 
1.12) are perturbed by the following geometric factors. 

1, ) The effects of changing ground resolution element size 
at different scan angles, and effectively recording on a 
cylindrical rather than plane surface, which causes a 
panoramic appearance. These are termed "scan angle effects"* 
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2. ) The effect of neglecting relief difference of object points 

above or below the assumed datum surface. This will be 
called "topographic effect". 

3. ) Removing the assumption of ideal flight conditions leads 

to "sensor exterior orientation effects". 

1*. ) Relaxing the assumption that each scan line is instantaneously 
recorded, but in fact, requires some finite recording time, 
leads to "scan time effects". 

l.U.2.1 Scan Angle Effects. Although each element along the 
scan in the image is displayed with equal width, significant dis- 
tortions occur along each scan due to a changing ground size of 
resolutions elements. From Figure l.l6, the element size in the 
direction of flight, or X direction, is given by 

e = 3(Z C - Zj) sec@ j ll.lU) 

J 

In the direction perpendicular to the direction of flight, 
or along-scan direction, the ground coordinate of the point J is 
given by 

Yj = Y c + (Z c - Zj) tan0 j (1.15) 

If the sampling angle, y, (see Figure l.l6) is sufficiently small, 
as is usually the case, then it can be considered as a differential 
change in the scan angle 0^, that is 


Y 


d0 


( 1 . 16 ) 
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Differentiating Equation (1.15) with respect to 0^ , the resulting 
expression is 

dY, * (Z c - Z 5 ) sec 2 0 1 d0 

or j J J U.iT) 

e Yj = ^ z c - se ° 2(1 j 

Figures 1.18 and 1.19 show plots of the ground resolution element 
size as a percentage of flying height above terrain in the X and Y 
directions, respectively. These plots show various values of the 
angular resolution. Most aircraft data handled by the LARS system 
is digitized in the y = 3-6 mrad range in the Y direction. The 
angular resolution of the ERIM M-7 scanner, from which most present 
LARS aircraft digital data is derived is on the order of 2-3 mrad, 
representing 3. The scan angle limits of the ERTS scanner is on 
the order of 5*75° either side of nadir. 

As an example consider a typical aircraft flight as utilized 
by LARS, with nominal flying height above terrain of 5000 ft. (1.52 km), 
digitization parameter y - 6 mrad, 3 = 2.5 mrad. Then the ground 
size of a resolution element would vary from 12.5 ft. x 30 ft. 

(3.81 m x 9*14 m) at nadir to 17.7 ft. x 60 ft. (5.39 m x 18.29 m) 
at a scan angle of 45°. As may be seen from Figures 1.18 and 1.19, 
the ground resolution element size increases very rapidly, particularly 
in the Y direction, for scan angles in excess of 45°, which would 
result in great image distortions. Therefore, most aircraft scanner 
data is collected within the range of a <t 40°. For small scan 
angles (e.g. ERTS, in which a < 6°) the effect is not as serious 
as for larger angles. 
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Figure 1.19. Resolution Element Ground Size in Y Direction 
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In addition to the distorting effect of changing ground reso- 
lution element size, there exists in each scan line a panoramic 
image displacement. This displacement is due to the fact that the 
data are effectively imaged upon a cylindrical surface, rather than 
a plane. From Figure 1.17b, this displacement may be seen to be 
the difference between some hypothetical image position, denoted a’, 
and the actual image position a. The resultant image displacement 
is given by 

6y = Oa' - Oa = c(tan6j - 6^) (l.l8) 

Figure 1.20 shows graphically the magnitude of this displacement 
relative to the effective principal distance, c, for varying scan 
angles. It should be noted that this displacement increases very 
rapidly for scan angles in excess of ^5°. 

l.U.2.2 Topographic Effect. The neglect of terrain variation 
is not uncommon in the interpretive treatment of multispectral 
scanner data arrays. The resulting displacement may become appreciable 
at large scan angles. Referring to Figure 1.17b, if point A is 
imaged and assumed to be on some datum surface , while in fact it 
lies at some elevation Zj above datum, the effect is to image the 
point as if its Y position was at A' on the datum. The effective 
ground displacement due to neglecting topography may be given by 

6Yj = Zj tan0j ( 1 . 19 ) 

in which Zj is the element elevation relative to datum, and 6Yj 
is the resulting Y displacement at ground scale. The proper image 
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position for this element would then be at a", if the element is 
to retain its proper Y position, after reduction to the datum surface. 
This error, 5Yj , will vary from zero at nadir to a value equal to 
Zj at a 45° scan angle and would increase rapidly for scan angles 
greater than 45°. In order to alleviate this error, some method 
must be chosen to assign elevations to every element with an acceptable 
accuracy. This problem will be dealt with in greater detail in 
a later section. At this point it should be mentioned, however, 
that if multiple coverage of ground areas are available, for example 
from side lapping flight lines, the possibility exists to obtain 
these element elevations directly from the imagery, by forming inter- 
sections. However, if only singly scanned data are available, then 
some source external to the data arrays must be used to provide 
these element heights. 

1.4. 2. 3 Exterior Orientation Effects. In addition to scan 
angle and topographic effects, a major source of error is that of 
sensor exterior orientation variations during scanning. In actuality, 
of course, an aircraft cannot be controlled such that its elevation 
and velocity remain constant, nor can it be constrained to lie 
on a perfectly straight flight path. In addition, an aircraft 
cannot be stabilized such that roll, pitch, and yaw are held negligible. 
Many scanning systems, however, are roll stabilized, such as the 
ERIM M-T scanner which supplies the present bulk of aircraft data 
to LARS. For these systems, the w term may be approximated by zero. 

Figures 1.21 and 1.22 illustrate the resulting deformations 
which occur when the perfect orientation assumptions of Section 1.4.1 
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are not realized. Figure 1.21 illustrates the effect of two of 
the angular orientation elements and one of the positional elements. 

In Figure 1.21a, the effect of pitch (<J>) variations down the strip 
is to cause the scanner to ’’gallop" in the X direction, alternately 
increasing and decreasing the nominal scan overlap value, S, described 
in Section 1.4.1. The effect of a change in flying altitude is 
to cause a spread of the sensed ground element in the direction 
normal to flight as may be seen from Figure 1.21b. 

The yaw (k) variation is due to the rotation of the scanner 
about a near vertical axis as it flies the strip. For example, 
in Figure 1.21c, assume a road is exactly perpendicular to the 
flight direction, but is imaged when the aircraft has a yaw. Then 
when the imagery is displayed, each line is displayed horizontally, 
and the x coordinates along the road will be displaced such that 
the road will no longer appear perpendicular to the flight line axis. 

Figure 1.22 [20] shows the resultant ground coverage of a scanner 
under various perturbations in sensor exterior orientation elements . 
From Figure 1.22, an important aspect of scanned imagery becomes 
apparent. Because only two-dimensional information is available 
from the imagery, it becomes impossible to separate the effects of 
each orientation element variation. For example, considering 1.22a, 
a constant aircraft pitch in scanning several successive lines is 
impossible to distinguish from the effect which would be recorded 
if a constant X shift of the coordinate system was done. From 
1.22i, it is not possible to distinguish between a linear rate of 
aircraft pitch, and a linear rate of change of aircraft velocity. 
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Figure 1.21. Effect of Some Aircraft Instabilities 
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Figure 1.22. Resultant Earth Coverage of Scanner in the Presence of 
Different Vehicle Orientation Irregularities (Adapted 
from [20]) 
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Thus, it would he unreasonable to expect the recovery of all orientation 
parameters from the imagery, especially if only singly scanned imagery 
is available. Kratky [21] also points out this inseparability of 
orientation elements based upon scanner imagery. 

The variation in these elements of sensor exterior orientation 
down the flight line are in general stochastic in nature, or random 
time functions. The cumulative effect of these, however, may be 
considered deterministic or analytical, particularly for short sections 
of the flight line [22]. It may be possible then to approximate 
these variations within a short section using such functions as 
polynomials or harmonics [21, 22, 23]. 

Ideally, the mathematical functions expressing the behavior 
of exterior orientation elements should be In terms of time. However, 
time may not be recorded with sufficient accuracy to be relied upon 
for such analysis. Alternatively, if the concept of a constant 
effective speed of "film" travel (v) is utilized, then the x coor- 
dinate of the imagery may replace time in these functions. When 
the data is in the form of digital arrays, x may further be replaced 
by the scan line number as the independent variable, 

x = v(T - T q ) (1.20) 

t = x (1.21) 

v 

where t = T - T Q . 
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As an example, if Y c is assumed to be approximated by a second 
order polynomial in time within a section of a flight line 

Y c = a’ 0 + a’j t + a' 2 t 2 (1.22) 

then the resulting polynomial form in terms of the image x coordinate 
would be 


Y c “ a Q + x + &2 x 2 


(1.23) 


in which 


a l = a V v 


(1.24) 


Derenyi and Konecny [22, 24, 25] investigated the effect of 
differential changes in each of the exterior orientation elements 
on the ground coordinates. The equations utilized were based upon 
Halle rt ’ s differential formulas [26] modified for scanned imagery. 
The resulting form of the equations is 


dXj “ (Z c - Zj) tanQj d< + (Z c - Zj)d<f> + dX c 
dYj = tan0j dZ c - (Z c - Zj)(l + tan 2 0j)du + dY c 


(1.25) 


For the present investigation, ground displacements were calculated 
for representative values of the orientation element differential 
changes for various combinations of flying height and scan angle. 
Figures 1.23, 1.24, 1.25, and 1.26 illustrate the ground displacements 
due to some of the terms in Equations 1.25 . Figure 1.23 depicts 
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Figure 1.23. Ground Displacement in X Direction 
Due to Change in Yaw (die) 
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Figure 1.2U. X Ground Displacement Due to Pitch Change (d$) 
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Figure 1.25. Y Ground Displacement Due to Flying Height Change (dZ c ) 


Scan Angle (0) 

Figure 1.26. Ground Displacement in Y Direction 
Due to Change in Roll (du>) 
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the X ground displacement relative to flying height above terrain 
due to a differential change in yaw (die). Figure 1.24 depicts the 
absolute displacement at ground scale for various values of flying 
height due to a small change in pitch ( d<J>) . Figure 1.25 illustrates 
the resulting Y coordinate change as a function of scan angle, 0, 
for various values of differential changes in flying height (dZ c ). 

Figure 1.26 gives the Y displacement relative to flying height as 
a function of scan angle, due to changes in the roll angle (dto). 

1.4.2. 4 Scan Time Effect. If the assumption that each scan 
line is recorded instantaneously is relaxed, then the data would 
be recorded as shown in Figure 1.27, assuming aircraft velocity, V, 
to be constant during the very short time span of a single scan. 

The Xj coordinate of any point may then be written as 

Xj = X c + yj t r V (1.26) 

2 ir c 

in which yj is the along-scan image distance to point j, or the 
column number for digitally recorded arrays . The t r term is , as 
before, the period of revolution of the scanning mirror, and c is 
the sensor constant (see Figures 1.16, 1.17). 

The period of revolution, t r , is a well defined term for most 
scanning systems. However, the aircraft velocity V, is a function 
of time, and some difficulty may be encountered in accurately monitoring 
this velocity. It is possible, however, to include this scan time 
effect into the k term of sensor exterior orientation as observed 
from Figure 1.27. 


Figure 1.27. Equivalent Yaw Due to Scan Time Effect 
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1.!? Review of Previous Investigations 

The imagery generated by an optical-mechanical scanner contains 
many of the features of imagery generated by more conventional strip 
and panoramic photography. Consequently, many of the geometric 
characteristics of the panoramic and strip imageries serve as a basis 
for analysis of scanned imageries, and the review will begin with 
a discussion of these more conventional sensing systems. This will 
be followed by a chronological review* of previously published investi- 
gations of the geometric aspects of non conventional sensors. These 
sensors include radar sensors such as Plan Position Indicator (PPI) 
and Side Looking Airborne Radar (SLAB), infrared line scanning systems 
(IRLS), and multi spectral scanning systems (MSS). The last portion 
of the review deals with the recent advent of remote sensing from 
spacecraft, and the geometric peculiarities of such systems, par- 
ticularly the simplifying assumptions which are possible due to the 
relatively high stability of such sensing platforms. Particular 
attention will be paid to the Earth Resources Technology Satellite 
(ERTS), which contains an MSS system. The inclusion of this sensing 
system led directly to several fruitful Investigations of scanner 
geometry. 

In the late 19th century an Italian named Porro developed a 
camera employing the panoramic principle. The instrument, fitted 
with a telescope and level, recorded photographs on a strip of 
sensitized paper held against the surface of an upright cylinder. 

Only fairly recently, however, have intensive efforts been directed 
toward uhe design, development, testing, and routine application 
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of panoramic cameras to aerial photography [2j]. Continuous strip 
photography was first used in 1932, and was initially intended for 
low altitude, high speed reconnaissance. 

Early papers "by Katz [28, 29] presented methods by which object 
heights could be extracted from continuous strip imagery using 
parallax measurements made on the photograph. Wohl and Stickle [ 30 ] 
in 1959 derived data concerning traffic velocities and volumes based 
upon the geometry of the continuous strip camera. In an article 
published in 1962, Elms [31] proposed the possibility of using the 
convergent strip camera system as a mapping tool. In that article, 
the requisite equations for determination of ground elevations and 
distances were derived under the assumption of the aircraft flying 
straight and level. The effects of errors due to insufficient 
image motion compensation, aircraft velocity and position variations, 
and aircraft angular orientation variations were investigated and 
their effect upon the original equations was presented. The conclusion 
which Elms reached was that mapping was possible using imagery from 
a strip camera, although the use of conventional plotting systems 
as they then existed was precluded. 

Ockert [32], in an article published in I960, compared the 
frame and strip cameras for potential use in satellites. The con- 
clusions he reached were that the frame camera, due to its highly 
refined and perfected geometry was preferable to the strip camera 
with its complex geometry. The problem of time dependent orientation 
parameters was recognized, with its attendant difficulties in performing 
radial and space triangulation based upon strip imagery. 
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In a paper published in 1961 [33] Abraham developed projection 
equations for frame, strip, and panoramic cameras, and used these 
equations to perform variance propagation investigations in an attempt 
to quantitatively evaluate the accuracies obtainable from the three 
camera types. After first formulating the projection equations 
of the frame camera, Abraham developed transformations for the pan- 
oramic and strip imageries to reduce measured coordinate values in 
these systems to equivalent frame camera coordinates . These trans- 
formations were introduced into the original projection equations 
to form projection equations for the panoramic and strip cameras. 
Variance propagation was carried out for a test case and numerical 
values computed. The results indicated that expected ground coordinate 
standard errors were about 10% greater for X coordinates obtained 
from panoramic images than those from frame photography, and about 
15% greater for Y coordinates. The strip camera yielded expected 
errors about lh% in X and 23 % in Y over those of the frame camera. 

A two article series [3*+» 35] authored by Itek Laboratories 
and published in Photogrammetric Engineering in 1961-62 did much 
to acquaint the photogrammetric community with the advantages, 
disadvantages and geometric peculiarities of the various panoramic 
cameras. The first article, published in Dec. 1961, dealt with the 
advantages of the panoramic camera (high resolution coupled with 
wide angle coverage), its disadvantages (longer cycling times resulting 
in unsuitability for low altitude, high speed aerial photography), 
and the two basic physical configurations possible in panoramic 
camera (direct scanning or rotating prism). The second of the 
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articles dealt with the geometric properties of the images produced 
by panoramic photography. In particular, effects of cylindrical 
imaging (previously discussed in Section 1.4. 2.1), the scan (sweep) 
positional effect due to camera motion during exposure, and the 
effect of image motion compensation were discussed. 

Kawachi [36] in 1965 presented image displacement equations 
for image motion due to angular orientation instabilities of frame 
and panoramic cameras. Kawachi also presented plots of rates of 
change for the roll, pitch and yaw orientation elements for a typical 
reconnaissance aircraft. These plots indicated the stochastic nature 
of the orientation elements . 

Gullicksen [37] in 1967 advocated the use of strip photography 
as a map substitute in order to eliminate the mosaicing necessary 
with frame photography. Equations for planimetric errors which 
resulted from the direct use of a strip photograph as a map substitute 
were presented, and a numerical example was computed. 

An article by Gill [ 38 ] in 1964 presented parallax equations 
for convergent panoramic photographs, and related these to a relative 
orientation procedure by analytic means. The photographs were 
relatively oriented in segments or parts, in order that the method 
could be used in an analytical plotter with limited format size. 

A test grid was simulated and the procedure tested using the AP-II 
analytical plotter, with the results indicating the feasibility 
of analytic relative orientation of panoramic photography. 

Hovey [39] in 1965* reiterated the advantages of panoramic 
photography, reviewed the image displacements inherent in the imagery. 
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and advocated, greater use of panoramic imagery by the photogrammetric 
community. 

Skiff [Uo] in 1967 published an article treating strip and 
panoramic photographs analytically. The article used tensor notation 
throughout. Collinearity equations were presented for frame, strip 
and panoramic photography, and the basic equations necessary for 
triangulation using two photos were obtained. Skiff included the 
concept of time varying orientation parameters, and assumed a linear 
time function for aircraft positional elements for the strip case, 
and a uniform angular velocity for the panoramic sweep. Derivatives 
of the resulting equations were presented for subsequent linearization. 
In concluding. Skiff suggested the use of a constant orientation 
assumption over a portion of the imagery. ' 

In a very important paper [4l] published in 19&7* Case formulated 
analytic expressions for projection and collinearity equations 
associated with panoramic and strip photography. The basic principle 
advocated in this paper is the reduction of image coordinates from 
strip or panoramic photography to coordinates on an "equivalent 
frame photograph", through simple transforming equations. The re- 
sulting image coordinates, it was believed, could then be treated 
using existing programs for space resection, relative orientation, 
and block adjustment. This formulation differed from that of the 
earlier work of Abraham [33] dealing with the same concept, in that 
a unified matrix approach was utilized, rather than deriving each 
case separately from the geometry of the projection. Derenyi [ 22 ] 
points out that Case advocated the use of a constant orientation 
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assumption for short sections of the imagery. This paper used 
standard matrix notation and a photogrammetric approach, and served 
as a basis for further analytical formulations of the less conventional 
sensors, such as SLAR and optical -mechanical scanners. 

A panoramic rectifier described by Wright [42] effectively 
projected a panoramic negative from a cylindrical surface having as 
its radius the principal distance of the imaging camera onto a tangent 
plane, forming a rectified image which would eliminate scan angle 
effects. Problems in lens design for this rectifier were presented 
in this paper. 

Masry [23] in 1969 published an article treating convergent 
stereo strip photography in a rigorous analytical manner. The neces- 
sary coplanarity condition for relative orientation was presented, 
and it was noted that all exterior orientation elements of the strip 
camera were dynamic in nature, that is, stochastic functions of 
time. Masry proceeded to make the assumption that the angular or- 
ientation elements of roll, pitch, and yaw were constant during any 
short period, or piece of imagery. The positional elements of 
relative orientation X c , Y c , Z Q were assumed to have some general 
polynomial form for both sides of the stereo strip photography, 
resulting in a polynomial form for the base elements of relative 
orientation B x , By, B z . Relative orientation was presented by 
dividing the strip into sections of a length approximately equal 
to the image length covered by the parallax angle. Corresponding 
sections were relatively oriented using a successive procedure in 
which the orientation elements of the section just relatively oriented, 
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on one side of the film, were used as the elements of orientation 
for the opposite section on the other side of the film, since they 
were simultaneously recorded. 

Masry then presented projection equations to form model coordinates 
if orientation parameters were assumed known, and test cases were 
run on both simulated and real imagery. Problems occurring were 
due primarily to the instability of the solution for orientation 
parameters and resulting model coordinates, caused mainly by small 
parallax angles in the convergent strip camera. Model coordinate 
errors were found tc be largest in the Y direction (perpendicular 
to flight direction) . Masry concluded that if monitoring of orientation 
elements were done , the model coordinate stability would naturally 
be greatly enhanced. 

In 1971 Derenyi [22] completed a thesis relating to relative 
orientation of strip imagery. The dynamic nature of the imagery 
caused by continuously changing orientation parameters was recognized 
and described. .Two basic schemes of relative orientation were pre- 
sented to cope with this problem. 

In the first of these, called the line-by-line method, two 
corresponding lines perpendicular to the flight direction were rel- 
atively oriented, and this procedure was repeated down the strip 
to form a series of line models. Intermediate values were obtained 
by interpolation by harmonics or polynomials. Since only three orien- 
tation parameters were recoverable by relatively orienting a line 
pair, an investigation was conducted into the effect of the omission 
of two orientation parameters upon the remaining three. The results. 


as indicated "by both theoretical error propagation studies and tests 
on real and simulated imagery, indicated that serious errors would 
be introduced unless the omitted parameters were monitored or sta- 
bilized to the same precision expected from the relative orientation 
itself. 

In the second method investigated, termed the section-by-section 
method, portions of the strip were treated as a unit for the relative 
orientation procedure. In this method, some assumptions concerning 
the behavior of the exterior orientation elements over the section 
are necessary. Two assumptions were investigated in detail. In 
the first of these, the constant orientation assumption, the aircraft 
orientation parameters were assumed to remain fixed for each section 
pair to be relatively oriented. Error propagation studies by Derenyi 
indicated this method could result in serious errors due to the low 
attainable precision of the d<j> term. Since Derenyi’ s analysis dealt 
primarily with convergent stereo systems, this low precision, coupled 
with a small parallax angle of convergence, would result in large 
expected errors of model point elevations formed by subsequent inter- 
sections. This problem could be alleviated by larger parallax 
angles, but if this is done, then longer sections are required, 
resulting in a loss of validity of the constant orientation assumption. 
An analysis was also carried out to assess the effect of this error, 
i.e. the assumption of constant orientation when, in fact, orientation 
is dynamic in nature. It was found that large errors in d(j> and die 
result from relatively small non-uniformities in orientation. 

A further manifestation of the constant orientation assumption within 



a section is also apparent. At section boundaries down the strip, 
discontinuities in orientation parameters, and hence, model coor- 
dinates, will appear. 

As an attempt to partially overcome these problems , a changing 
orientation assumption, to be approximated by linear orientation 
changes within a section, was also investigated. It was found that 
significant improvement in accuracies resulted from this more realistic 
assumption, and the discontinuities at section born iaries disappeared. 

As a further attempt to strengthen relative orientation procedures, 
Derenyi suggested a "triple channel recording" system, in which 
a convergent stereoscopic system would include not only fore and 
aft looking sensors, but also an additional vertical one. By this 
method, the three lines recorded could be treated much as lines 
from a single frame in a frame camera, and many of the techniques 
of frame photography could be employed. It was felt that use of 
this scheme would allow greater parallax angles, and hence help 
alleviate the low precision attainable in the d<i> term, and subsequent 
model point elevations. 

Derenyi also analyzed the sidelapping flight line configuration. 

As mentioned in Section 1.4. 2. 3, the inability to separate analytically 
between some orientation elements results in the possibility of 
recovering only two orientation elements, and then only by using the 
line-by-line approach. Derenyi draws the conclusion that this flight 
configuration was therefore unsuitable for relative orientation. 

It should be mentioned, however, that this means relative orientation 
in the strict photogrammetric sense, in which it is desired to recover 



five orientation elements. While this is clearly impossible using 
sidelapping flight lines, it may be possible, after some type of 
relative or absolute orientation has been done, to obtain model 
elevations by intersection. 

In fact, Leberl [6] in 1972 studied the propagation of error 
from orientation parameters to model coordinates for optical-mechanical 
scanners. His results indicated that, with regard to error propagation, 
the sidelapping of flight lines was the optimum scheme in order to 
minimize cofactors of the Z model coordinate , or elevations . 

A recent use of panoramic photography of general interest was 
reported by Alderman [43], in which a convergent panoramic photography 
system was used in establishing control on the lunar surface . Oper- 
ating from the Apollo 15 Scientific Instrument Module, an optical 
bar camera was used to generate overlapping panoramic photography. 

The resulting imagery was used to generate models by the process 
of relative orientation utilizing the constraints imposed by the 
orbit of the spacecraft. These were subsequently assembled into 
a strip using the method of Schut [44]. Subsequent transformation 
of the assembled strip to control points was then carried out. 

A linear rate of change for exterior orientation elements was assumed 
to accomodate the dynamic nature of the imagery. 

With the advent of the concept of "remote sensing" in the general 
sense, beginning in the early 196q's, the photogrammetric community 
has been faced with problems of an unusual and challenging nature. 

The unconventional sensor types advocated for use , such as SLAR 
and MSB, although having unique data gathering capabilities such 


as all weather operation and image forming capabilities outside 
the visible range of the spectrum, suffer from problems of reduced 
resolution and image distortion due to the dynamic mode of image 
formation. 

The use of radar as a remote sensing system has been progressing 
for some time, beginning with the Plan Position Indicator (PPI) 
first developed around 1940, with Side Looking Airborne Radar (SLAB) 
coming into existance around the mid 1950's. Leberl [45] presents 
a detailed history and literature review on radar remote sensing 
systems . 

Suits [46] reported in i960 on the nature of infrared radiation 
as a source of imagery. Recording systems based upon both cameras 
and scanner-detector configurations were presented. Colwell [47] 
in 1963 published an article which did much to acquaint the photo- 
grammetric community with the basic matter-energy relationships 
necessary for image formation by remote sensing. Leonardo [48] 
in 1964 published a further article dealing with general advantages 
and problems associated with remote sensing systems. 

In an early paper, Harris and Woodbridge [49] presented the 
basic configuration for a thermal scanning system and the concepts 
of spectral radiance in the emissive region of the spectrum, as 
well as atmospheric transmission. 

Suits [50] in 1966 advocated that security restrictions on 
infrared scanning devices be lifted, in order that the scientific 
community as a whole might benefit from this relatively new data 
gathering system. Suits reported on an anticipated reduction in 
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the classification of instruments, thus alerting the photogrammetric 
profession to its availability. 

Other sources of general information on remote sensing and 
remote sensing systems may be found in references [4, 10, 11, 27, 

51, 52], as well as the series of LARS annual reports [53-56], 

and the proceedings of symposia held at the University of Michigan [57]- 

In the mid 1960's, the airborne multispectral scanner was used 
at the Willow Run Laboratory of the University of Michigan (now 
Environmental Research Institute of Michigan). Since 1966 Purdue 
University's Laboratory for Applications of Remote Sensing has been 
formulating and' developing a digital system of MSS data handling to 
perforin automated interpretation based upon spectral differentiation 
by statistical and pattern recognition techniques [53, 54, 55, 56 ] . 
Recent papers by Wilson [58] and Lapides [59] have emphasized the 
utility of multispectral scanners and enumerated their unique advan- 
tages, as well as limitations associated with the data available 
from them. Carnes [6o] reported on advances in auxiliarly equipment 
associated with such scanning systems. 

Until recently, much effort has been spent on investigations 
into the radiometric and interpretative aspects of data gathering, 
data handling, and automated interpretation. The geometric aspects 
of the data arrays generated were ignored or sublimated to the more 
pressing problems of research into and development of the radiometric 
and interpretive aspects noted above. As the analysis of non-metric 
facets progressed, however, the problems of geometric factors in these 
scanning systems became more apparent and woric in this area was begun. 
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In a series of two papers [24, 25] Derenyi and Konecny first 
applied a photogrammetric approach to the optical-mechanical scanning 
system. These papers dealt with the basic recording system, the 
scan angle effects of varying resolution size and panoramic displace- 
ment mentioned in Section 1.4. 2.1, and the dynamic nature of scale 
factors involved in a scanning system. The study of exterior orien- 
tation effects was based upon a modification of Hallert's differential 
formulas. The effect of errors in determination of angular resolution 
size was also assessed. In the second of the papers, the authors 
presented a projection equation in matrix form which considered the 
dynamic nature of the imagery, accounting for the movement of the 
projection center during recording. 

Konecny, in a series of papers, [7» 6l, 6 2, 63 ] summarized the 
basic projective relationships associated with various remote sensing 
systems. These papers dealt with methods of approximate relative 
orientation and rectification of imageries recorded with dynamic 
systems and discussed accuracies which may be expected utilizing 
these systems. The latest paper [ 7 ] provides an excellent summary 
of most of the work done on geometry of non “Conventional remote 
sensors to date. 

A paper by Taylor [64] provides useful insight into rectification 
procedures for infrared line scanning sensors. Rectification equations 
are derived from basic projection equations by first assuming a 
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perfect orientation case, and analyzing image displacements after 
subsequent perturbations in roll ( to) , pitch ($), and yaw (k). 
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A series of papers [6, 45, 65, 66, 67, 68] written "by Leberl 
treat geometric aspects of the two most popular non-conventional 
remote sensors, side looking airborne radar (SLAB) and the optical- 
mechanical scanner. Leberl was primarily interested in SLAB imagery, 
and his work culminated in a doctoral thesis [45] which is particularly 
notable for two reasons. The first of these is the extensive bibli- 
ography included on radar in general, and in particular, SLAR. 

The second notable point in this work is the introduction of inter- 
polative techniques for the rectification or restitution of SLAR 
imagery. Using stochastic concepts, these interpolative techniques 
present a viable alternative to the projection equations for image 
rectification or restitution. These techniques may prove to be 
quite useful for the rectification of digital MSS data arrays, since 
they appear well suited to this digital data handling format (parti- 
cularly when scanned terrain is relatively flat). 

A paper by Markarian, et. al. [69] in 1971 presented a scheme 
of utilizing general polynomials for rectification of digitized 
imagery. The concept of treating the data as an uncorrected input 
array used to generate a rectified output array was utilized, and 
the principal of digitally "stepping" on the output array and placing 
gray scale values from computed positions in the distorted input 
array was advocated. Further, it was recommended that the "nearest 
neighbor" be used as the gray scale value from the input array, 
rather than to attempt interpolation of gray scale values. 

Increasing interest in these geometric aspects of remote sensing 
systems is evidenced by the increase in literature in this area 
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very recently. Masry and Gibbons [TO] reported on the use of the 
analytical plotter for rectification of scanned infrared imagery. 
Rectification of planimetry using the LR-2 line drawing rectifier, 
as reported by Forrest [71] handles either SLAR or scanned imagery. 

In a paper by Bosman, et. al., [72] a programming system called 
KARIN is described for planimetric mapping from single or overlapping 
strips of remote sensing data such as SLAR or infrared line scanner 
(IRLS). Graham [73] and Yoritomo [74] described restitution schemes 
for radar imagery at this same meeting. 

In the early 1970’s the use of satellites and spacecraft as 
remote sensing platforms was realized. The culmination of this long 
awaited event came about with the launching of the ERTS-1 satellite 
on July 23, 1972, which contained both return beam vidicon (HBV) 
and multispectral scanner (MSS) imaging systems. Because of this 
event, considerable effort was expended prior to, and immediately 
after, the launch on geometric problems associated with the MSS. 

The Skylab manned orbiting Earth resources laboratory contains 
photographic and MSS imaging systems . The Skylab MSS is of the 
conical scanning type. 

General information on the imaging of the Earth’s surface from 
space , and the ERTS-1 and Skylab systems and their sensors in par- 
ticular may be found in references [75, 76, 77, 78, 79, 80, 81, 

82, 83 ]. In a very revealing paper, Forrest [75] explains the unique 
advantages of using satellite imagery, particularly from scanners, 
for cartographic use. Consideration was given to the scanner geometry 
and resolution, and the use of ground control points extracted from 
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existing maps was advocated to resect for sensor orientation elements 
for the precision required for mapping beyond that obtainable from 
satellite ephemeris data. Forrest points out the greater stability 
of orientation in spacecraft over that of airborne sensors. The 
conclusion was drawn that mapping from space images economically 
compares very favorably with that of conventional techniques, and 
this cost may be spread among users of the imagery outside of those 
with cartographic interests. 

In two recent papers [8l, 82], Doyle summarizes the types of 
sensors which have to date been used from spacecraft, and describes 
some of their uses. Processing of the data and imagery generated 
is also discussed. 

The geometric problems associated with imaging with spacecraft 
mounted sensors have been addressed by several investigators. 
Although, as previously mentioned, spacecraft orientation insta- 
bilities are generally less pronounced in these systems than with 
airborne systems, these investigations yield valuable insight into 
possible analysis methods for airborne systems. 

Colvocoresses [84] as early as 1970 published an article dealing 
with geometric errors to be expected in ERTS imagery. He analyzed 
image displacements for 5 possible sources: Earth curvature, at- 

mospheric refraction, camera obliquity (or exterior orientation 
angular elements), terrain relief, and map projection error. The 
effect of atmospheric refraction was found to be negligible , and 
that due to terrain relief would be minimal except in extreme cases. 
Other errors were easily modeled except for the obliquity error. 
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Colvocoresses advocated not processing images in an element by element 
manner except for special purposes. 

A series of papers by Kratky [21, 85, 86, 87] addressed these 
problems from the photograrametric viewpoint. Of particular interest 
for this investigation is the formulation of functional forms to 
account for orientation instabilities of the scanning system. 

After presenting the basic projection and collinearity equations 
for 'ihe scanned imagery, the use of two alternate functional forms 
to approximate orientation variations with time is investigated. 

These are polynomial and harmonic functions. Kratky presents the 
equations in matrix form for resection to solve for the parameters 
of these functions, including explicit presentation of the coefficient 
matrix and constant vector for each of the functional forms considered, 
when substituted into the linearized collinearity equations. Because 
of the smaller and smoother orientation behavior of spacecraft systems 
over airborne systems, it was possible to incorporate small angle 
approximations into the analysis, yielding simplified equation forms. 
Kratky also poiL-t’t out the inability to distinguish between the 
correlated orientation elements and X c , the pitch and sensor X 
position, using singly scanned imagery. Thus, it is indicated, 
resection for all orientation parameters is impossible. However, 
the remaining parameters may compensate for these neglected elements, 
and result in acceptable image positions after subsequent projection. 

The geometric processing of ERTS imagery in the United States 
has been discussed by several investigators [20, 88, 89 ]. Forrest [20] 
describes the bulk image processing scheme (currently called -"system 
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corrections"), as well as the more refined precision processing 
("scene corrections") used on selected imagery as requested "by the 
user. The bulk image processor transforms the scanned imagery to 
a central perspective projection. The purpose of this transformation 
is to generate a series of "frames" from sections of the imagery 
as they would appear if recorded by a frame camera. The element 
corrections are fully calculated for intersections of a nine by nine 
grid. These displacement values calculated incorporate corrections 
for Earth rotation, as well as roll, pitch, and yaw values, which 
are taken from satellite data. Image positions between these 9 by 
9 array points are assigned by linear interpolation. 

The precision processed (scene corrected) images are generated 
as a subset of the above bulk processed (system corrected) imagery, 
incorporating more sophisticated techniques for geometric improvement. 
The use of ground control points extracted from existing maps, and 
identifiable on the imagery are used to refine image position cor- 
rections. Forrest points out several advantages in the use of ground 
control points. 

,1. ) The use of control points allows for obtaining the positional 
accuracy desired in the precision processed imagery. 

2.) Failure of the attitude sensor (the sensing unit which 

monitors the angular orientation elements of the spacecraft 
sensor) will not jeopardize the possibility of geometric 
improvement of the imagery. 
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3.) The direct determination of sensor orientation elements 
by resection from these points serves as a check on the 
attitude sensor, and serves to refine the orientation element 
values used in the bulk processing. 
h.) The use of redundant control points and least squares 
adjustment techniques allow the estimation of obtainable 
accuracies for orientation elements and indicate the positional 
accuracies of the control points used. 

The resection equations used to correct the MSS imagery in 
precision processing solve for eight parameters. These parameters 
reflect the image displacements caused primarily by errors in pitch, 
roll, and yaw, and by rates of change of these orientation elements. 

Once again, due to the small, smooth nature of the spacecraft orien- 
tation variations, the use of polynomial functional forms may be 
used to represent these image displacements . The equations used are 

Ax = a Q + [l + (x/H ) 2 ] + a^x + a^y + a-xy 

(1.27) 

Ay = b Q + [1 + (y/H s ) 2 ] + b-jx + b 2 y + b^xy 

in which Ax and Ay represent image shifts from the x and y positions 
recorded on the bulk processed imagery, H s is the spacecraft altitude 
at image scale, and a^, b^ are the eight parameters to be solved 
for in the spatial resection. It should be noted that terrain 
relief is neglected in these equations. ICratky [21] points out that 
the effect of terrain relief is less than the ground resolution 
size of an element unless terrain elevation differences on the order 
of 2600 ft. (800 m) are encountered within a single frame. Such 
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relief differences do exist in some ERTS frames, and the planimetric 
accuracy in such frames will suffer accordingly. 

Since digital data handling techniques for MSS data have become 
operational at LARS and elsewhere, it is desirable to investigate 
the geometry of such imagery in greater detail. The geometry of the 
scanning system will be analyzed in this investigation, with particular 
emphasis to be placed upon techniques which are suited to the digital 
data handling approach. Existing published literature in this area 
has dealt with data from spacecraft platforms, in which simplifying 
assumptions may be made, or for digitized frame photography. The 
investigation in this thesis will attempt to deal with the more 
general problem of gathering MSS data from an airborne platform, 
for which the simplifying assumptions based upon platform stability 
do not in general hold. This combination of a general geometric 
analysis, valid for both airborne and spaceborne systems, coupled 
with the digital data analysis, it is felt, will provide a much 
needed work not present in existing literature. The chapter immediately 
succeeding this will therefore deal with the general geometric 
problems associated with the MSS scanning system, in order that the 
relationships obtained may be utilized for the analysis of specific 
areas of interest in MSS digital data systems , 
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2. GEOMETRIC ANALYSIS OF MSS DATA ARRAYS FROM SINGLY SCANNED AREAS 

2.1 Remote Sensing as a Transformation 
According to Mikhail and Baker [5] s remote sensing may be thought 
of as a mapping of multidimensional object space onto another space, 
the sensor space , having the same or fever dimensions . Such a mapping 
may be effected through a transformation taking a multidimensional 
object space vector into multidimensional sensor space vector to 
be stored for subsequent retrieval and analysis. This may be termed 
the ’’remote sensing transformation”. The dimensions of each vector 
vill depend upon the degree of simplification of the physical phenomena 
involved, and the characteristics of the sensor employed. Examples 
of information content which may be included in the object space 
vector are [ 5 ]: spectral radiance from each resolution element in 

a particular wavelength band, polarization of the radiant energy with 
respect to the object space coordinate system, and coherence, both 
spatial and temporal, of the radiant waves. Components of the sensor 
space vector may include: spectral irradiance incident on the sensor 

in each wavelength band considered, the direction of polarization 
(if retained) with respect to the sensor coordinate system, the 
degree to which coherence is attenuated, and frequency and phase 
shifts of the carrier or force field if active sensing systems are 


utilized. 
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An important characteristic of sensors with image forming 
capabilities is that geometric image positions can be generated from 
the sensor space vector record stored. In most non-mapping systems 
to date , this geometric image position is not considered the primary 
goal of the recording. However, the fact that it can be recovered 
affords a unique opportunity to the interpretive data analyst and 
photogrammetrist. This leads to a unified data analysis concept 
in which analysis of both the non-metric data for interpretation 
as well as the metric record inherent in the data becomes possible. 
As indicated in Chapter 1 in the review of previous investigations , 
a two step sequential procedure has been used for some scanning 
systems (notably spacecraft systems) in which scanned data has been 
processed for geometric restitution, with the resulting data then 
utilized for automated interpretation. However, the possibility of 
simultaneous digital analysis of both the interpretive and geometric 
aspects has not been greatly exploited to date. With the more 
traditional photographic sensors, the nature of the methods used for 
both interpretation and photogramme try has largely precluded this 
unifying concept. The photo interpreter has used methods of visual 
inspection, relying on visual "keys" to arrive at decisions. The 
photogrammetrist has relied heavily upon accurate reconstruction of 
the geometry of the photograph to extract positional and metric 
information about the object space. An examination of the object 
and sensor space vectors for black and white frame photography may 
perhaps explain why these traditional fields of interpretation and 
photogrammetry remain separated. For each point imaged in the 
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photograph, the object space and sensor space vectors may be thought 
of as each having two parts. The first part will contain the non-metric 
information which may be sensed and recorded (such as spectral 
information, polarization, etc.), the second portion will contain 
geometric information. Consequently, the object space vector for 
each point may take the form: 


S 1 


i 



X 

Y 

Z 


spectral radiance from 
n wavelength bands within 
the photographic region 


point position in object 
space 


where consideration is given only to spectral values for the non- 
metric portion. The resulting sensor space vector is of the form 

film density recorded 


D 


for the point 


x point position in 

y image space 

These vectors are related through some physical-geometric transfor- 
mation which transforms the (n + 3) dimensional vector for a point 
in object space into a 3-dimensional vector in the sensor space. 

The physical transformation takes the n spectral values emanating 
from the point into a single density value D. The geometric portion 
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of the transformation takes the object space position (X, Y, Z) 
into the two-dimensional image position (x, y). The parameters 
associated with this geometric transformation are well known in 
photograrametry , and consist of sensor position (X c , Y c , Z ^) , sensor 
attitude (cu, a <)» and internal sensor characteristics (x Q , y Q , f, 
particularly for frame photography) . 

The problem with this transformation in general, which led to 
the traditional separation of the fields of photo interpretation 
and photograrametry, lies in the suitability of each part of the 
transformation for numeric analysis. Traditionally, the interpreter 
does not look only at a single point density value, but notes the 
density, texture, size, shape, shadow, tone, color, and pattern of 
objects before concluding their nature [90]. Recent attempts at 
automating this procedure, in which photographic densities are 
quantified and computations performed upon them are actually attempts 
to simulate these visual "keys". 

The geometric portion of the transformation, however, lends 
itself quite well to both analog as well as analytical treatments. 

This characteristic has spawned the practices of photo grammetry. 

However, with the advent of the unconventional sensors, par- 
ticularly the multispectral scanners, the extensive volume of the 
data severely curtails the role of direct human interpretation. 
Furthermore, the availability of the acquired data in digital form 
makes it possible to apply efficient computational decision making 
techniques and automated interpretation, as briefly described in 
Section 1.3. For the MSS system, considering only spectral information 
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for the non-metric portions of the vectors, the total object space 
vector would be of the form 

spectral radiance reflected 
and emmitted in n wavelength 
bands 

element position in object 
space 

time 

and the form of the resulting sensor space vector would be 

spectral irradiance recorded 
in m channels (m < n) 

element position in image 
space 

For the MSS transformation, the nature of the sensor makes it 
necessary to regard time as part of the object space vector, and the 
dynamic nature of the sensing process is reflected in the way in 
which the transformation parameters of the geometric portion are 
handled. 

Becaust? the sensor is recording continuously during the aircraft 
flight, the sensor position (X c , Y c , Z c )as well as attitude (w, k) 


f S’i 


! 



X 


y 


Si 



X 

Y 

Z 

T 



92 


become continuous stochastic functions of time. Other parameters 
of transformation for the geometric portion are c, a sensor internal 
constant analagous to the principal distance of panoramic recordings 
and v, a rate of imaging recording surface advance (see Section 1.U.2). 
Because of the dynamic nature of the sensing process , geometric anal- 
ysis becomes more involved and less precise. However, the main 
advantage of the system is the retention of the several spectral 
irradiance numbers for each element, which allows measurement and 
computation for interpretive purposes. 

In summary then, the MSS sensor allows the retention of sufficient 
information to make feasible a simultaneous analysis of both the 
physical (non-metric) characteristics as well as the metric (geometric) 
aspects of the data recorded. This simultaneous analysis may be 
brought about by combining the automated classification techniques 
discussed briefly in Section 1.3 with the geometric analysis methods 
to be discussed in this thesis. Therefore, the opportunity is avail- 
able for considering a unification of the activities of both inter- 
pretation and photogrammetry. 

This has not as yet been effected, however, for two primary 
reasons : 

1.) Historically, interpretation and photogrammetry have been 
considered separately. It was thus only natural to continue 
this trend even after these new sensors became available. 

Also, there exists a shortage of personnel who are adequately 
trained in both interpretive (particularly automated inter- 
pretation) and photogrammetric techniques. 
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2.) The concept of combining the two computationally oriented 
fields of automated interpretation and photogramme try, while 
logical, is quite a complex problem. Perhaps correctly, 
researchers to date have concentrated on one or the other 
aspect, in order to simplify the problems confronted and 

U 

gain new knowlege. 

The work done in automated interpretation has largely ignored the 
geometric aspects of the data vectors. As more sophisticated classi- 
fication algorithms became available, however, and as progress 
was made in automated interpretation , it has become apparent that 
the geometric aspects are gaining in importance. It is now becoming 
the task of the data analyst to raise not only the question of 
"what" (interpretation) but also of "where" and "how much" (geometry). 

This problem is well illustrated by the operational ERTS system, 
in which it was felt that, even for bulk (system corrected) images, 
some geometric analysis and preprocessing was necessary to efficiently 
utilize the data. It is therefore the objective of the work reported 
on herein to analyze in detail the geometric aspects of MSS digital 

t: 

data arrays in an attempt to provide answers to these queries. 

The subsequent section will introduce the basic transformation 
relationships which serve as the basis of analysis of MSS data from 
the geometric standpoint. 





9b 

2.2 T he MSS Geometric Transformation 
Figure 2.1 represents the recording of one channel of data 
during the scanning process of the MBS, in which the data array 
is considered as an image. Each array position would then have an 
associated "picture element" or "pixel". The two dimensional matrix 
of numbers would then be equivalent to the size of the imagery record. 
The line number i would represent a single scan, and within each 
scan the value j would represent the element position. With these 
comments in mind. Figure 1.17 may be thought of as the recording 
of the j-th element along the i-th scan line. Each picture element 
is recorded at a particular instant of time, thus making time an 
important variable. Following are the basic variables for subsequent 
analysis. Many of these variables have been previously introduced, 
but are included here for ease of reference. 

T q the epoch, or time instant, of beginning of 

recording (at zero x*-coordinatej 
t r time period for one revolution of the scanner 

T^j time instant of recording point 3 on scan 

line i 

0j instantaneous scan angle with respect to 

a scanning axis passing through the effective 
perspective center and the scan line center 
2 a total scan angle (a on either side of scanning 

axis) 

c constant representing a "principal distance" 

or the equivalent of the radius of a cylindrical 


recording drum 



Scan Line Number 
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v 


(x*. y*) i j 


y J 

Mij(t) 


k V kj 


(x c , x 0 , Z c ) 
(X, Y, Z)j 


speed representing ’'film" advance in case 
of pictorial recording, or its equivalent when 
recording in other modes (magnetic tape, 
etc.),', taken as a constant 
cartesian coordinates of point (either 
image or its equivalent in digital recording) 
as shown in Figure 2.1 

distance along scan line i from scan line 
center to 

instant aneous orientation matrix at 
(taking object system to sensor system) 
distance along the image center line from the 
beginning of the recording, represented in 
digital recording by the scan line number i 
scale factors between object and sensor spaces 
for point j on scan line i, (k^ = k ? j/cos Oj). 
k'j is the actual scale factor of panoramic 
recording. Referring to Figure 1.1Tb, this 
scale factor would be represented by the 
ratio of lengths Ca and CA. ,kj is the scale 
factor for the equivalent frame photograph, 
instantaneous position of exposure station 
at T i( j 

object coordinates of any point j on scan 


line i 
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Using the above variables, several useful relationships may be derived. 
If, as is usually the case, the period of revolution of the scanner 
(t r ) is a well defined quantity, then the time of imaging each element 
within the array may be determined from its array position by the 
relationship 



T 0 + (i-l)t r + fq + 6j) 


( 2 . 1 ) 


in which the substitution = 

The x* cartesian coordinate of the point may then be determined 


y^/c has been made. 


by introducing v defined above: 



The y* cartesian coordinate may be determined using the Pythagorean 
Theorem from Figure 2.1. 



It may be seen from Figure 1.17b that in the y direction, 
imaging takes place from about a perspective center having object 
space coordinates (X c> Y c , Z c ) at the instant of imaging point Pj. 
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The vector between the perspective center and a point in the object 
space coordinate system may then be written 


4 

A 



(2.5) 


An image space coordinate system may be defined by a right hand 
system in which the xy plane is tangent to the effective panoramic 
cylinder at the scan line center with the y axis in the direction 
of a scan line. The z image axis would be perpendicular to the xy 
plane and complete the right hand system. The vector above may then 
be written in this coordinate system as 


A = 


= Hu<t) 


X J - X c 
Y J - Y c 
Z J ' Z c 


( 2 . 6 ) 


In the image system, the vector from the perspective center to 
the image position may be written 


0 

c sin 0 
— c cos 0 


(2.7) 


The vector A would be represented in Figure 1.17b as the li-ne segment 


CA, and vector a would be represented by line segment Ca. 
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Assuming no systematic effects , these vectors must be coll inear . 

If vector a of (2.7) is multiplied by the proper scale factor, the 
resulting vector will be equal to that from Equation (2.6), leading 
to: 



0 


I 

x° >h 

1 1 

X >-t 

! 

1 

c sin e. 

= M. , (t) 

k’ . 

J 

-ij 

J c 

J 

-C COS 0 j 


Z i ■“!& 

J c 


( 2 . 8 ) 


Combination of the expressions of (2.2) and (2.8) with k'.^ = cos 0j 
results in the following four dimensional metric MSS transformation 
for rectilinear scanning: 


X 


o 

E-t 

1 

•H 

0 


v 0 


X - X 


= 

1,"3 


J c 

c/kj tan 


0 Si 1 ” 


- T = 

- c/k j 


- _ 


Z J - Z c 


This general transformation may then be used in subsequent analyses 
to investigate the geometric characteristics of digital MSS data 
arrays . 

When digital arrays are displayed, all elements in a single 
scan line i are displayed with a constant strip coordinate x., . 

The data may then be analyzed as if all points in a single scan 
line were imaged instantaneously with some residual yaw, as explained 
in Section I.4.2.U. The instant of recording would be assumed equal 
for all points j in a scan line i, i.e. x*^ = and T.. j = in 
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the first equation of (2.9)* The image coordinates which would be 
measured from a digital array display would be the i, j array element 
positions, or x^, y^ if the array is considered as an image. 

It may be shown that the projection equations presented in 
Section 1.4.1, which were generated from the basic geometry of scan- 
ning under ideal conditions, may also be derived from this general 
transformation if appropriate assumptions are made. For example, 
the first equation of (2.9) yields directly the transformation from 
time to image x coordinate as given previously in Equation (1.20) 
Under the assumption of ideal orientation, the matrix M..(t) 
becomes identity. If this substitution is made, the second equation 
of (2.9) yields directly the relationship of Equation (l.ll). 

Also under this assumption, if the last equation is solved for 
-c/k^j , and the result is substituted into the third equation, the 
resulting expression will be identical to Equation (1.12). 

The general transformation of Equation (2.9) may also be used 
to investigate the effect of small perturbations in orientation 
elements from the idealized case. Considering a single scan line, 
using = y^/c 821(1 introducing 


h j = c/k j 


( 2 . 10 ) 


the three dimensional geometric MSS transformation, after inversion, 


may be written as 
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“1 


hj tan 0j 


-h 




and rearranging, 


V 

Y, 


"x 

c 

Y_ 

+ 

0 

h. tan 0. 

J 


c 

-V 


i i 

- - b i - 


( 2 . 11 ) 


( 2 . 12 ) 


In order to investigate the effects of small changes in the 
orientation parameters X c , Y , Z c , w, <j>, k, on the coordinates of 
an object point the total differential of (2.12) must he evaluated. 

In order to evaluate the total differential, the partial derivatives 
of the equation set must he evaluated for each of the six orientation 
elements. Considering the sequential set of angles (w, <J>, k), 
according to Lucas [91] 



3(f) 


M. 



0 

0 

-1 



sm w 


-COS U> 



-sin w 
cos w 


0 

0 


0 

0 
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0 0 


0 0 0 


It is also necessary to evaluate the differential of the equations 
with respect to the parameter h^ , as it contains the scale factor 
kj which is a variable for each point. Proceeding with the evaluation 
of this total differential: 


+ M j tan ah. 


0 0 


0 0-1 


M j hj tan 0^ dco 


-sin (i) 


sin w 0 


-cos u> 0 


. j h. tan 0, d<}> 


+ M 


0-10 


0 0 0 


hj tan 0j die 


It is now assumed that the changes in angular orientation elements 


take place from a perfect orientation position of w = = k - 0. 


Then = I_, for which hj = Z c - Zj as may "be seen from the last 
equation of (2.11) under this ideal orientation assumption, and the 
resulting equations become 


~dx~ 


~dX~* 


0 


0 

dY 

= 

dY 

+ 

tan 0^ 

dh. + 

j 

h 5 

_dZ_ > 

j 

_dZ_ 

c 

_ -1 


h. tan 6. 
j d ■ 


(2.15) 



'-V 


-hj tan ©j - ] 

+ 

l 

0 o 

1 

dcjj + 

1 

o o 

1 


It should be noted here that Equation (2.15) represents a projection 
from two-dimensional space in the sensor to three dimensional object 
space. It is not possible from single imagery to determine all three 
coordinates (X, Y, Z)^, because the unknown scale factor k^ is 
implicit in h . Therefore, in order to derive the two possible 
coordinates it will be assumed that Z is constant for all j . Since 
Z is assumed constant , dZ ^ = 0 , and the last of the equations in 
(2.15) may be written explicitly as 


0 = dZ„ - dh. + hi tan 0. dto 
c j J j 


Then 
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Substituting this expression into the first and second equations 
of (2.15) results in 

dXx = dX_ - h, d<f> - h. tan 6, dx 

3 3 J 3 , (2.16) 

dYj = dY c + tan Gj dZ c + hj(l + tan Gj) dm 

Equations (2.l6) may he seen to he termwise identical to those 
of Equations (l,,25) in which h. - Z - Z. The contradiction in 

j c 

signs for the dw, d<J> , and dx terms is due to the fact that Equations 
u.25) were written as a special case of Hallert's original projection 
equations. Hallert, in his derivation [26], defined the positive 
sense of rotation for the angular elements in the opposite sense 
from that customarily used in the United States, and applied here. 

2.3 Collinearity Equations for MSS 
The collinearity condition introduced in Section 2.2 serves 
as the basis for collinearity equations. The collinearity equations 
for frame photography have been used extensively in the past as the 
basis of analytical photogrammetry [93, 9^1* Similar equations 
for scanned imagery may be derived using Equations (2.9). Considering 
a single picture element the i, j subscripts may be dropped and the 
geometric portion of Equations (2.9) may be written as 


0 


ffi ll m 12 m 13 


"x - x c - 

c/k tan 6 

= 

“21 ”22 “23 


Y -T 0 

_ -c/k _ 


J°31 m 32 m 33- 


- z - V 


(2.17) 
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Performing the matrix multiplication, dividing the first two 
equations by the third, and rearranging yields 


0 = 


m-^2_(X — X^) + m^g(Y — Y^,) + m-^^CZ — Z^) 


0 = tan© + 


(X - X c ) + (Y - Y c ) + m 33 (z - Z c ) 
m 2l (X - X c ) + niggCY - Y c ) + m^Z 


- Z„) 


(2.18) 


m 


3 l(x - X c ) + m 32(Y" Y c ) + m 33 (Z - Z c ) 


An alternate form may be derived using the transpose of (2. IT), 


0 = 


0 = 


(x c 


- X) + (Z - Z c ) tan 6 - 11132. 


0123 tan G - 11133 


-Y)+(Z-Z C ) mQ 2 tan 6 - 


(2.19) 


m. 


L 32 


If, in these equations, an ideal orientation assumption is made 
in which M is taken equal to I_ then Equations (2.19) reduce to the 
same results given in Equations (l.ll) and (1.12) in Section l.H.l. 

If six orientation elements were assumed unknown for each 
scan, it would be impossible to have a solution and derive metric 
information from MSS data arrays. Therefore, some type of functional 
behavior (polynomials, harmonics, etc.) must be assumed for those 
orientation elements. Once decided upon, object space control may 
be used to determine the coefficients of these functions (resection), 
then Equations (2.19) may be used to determine other object points 
(intersection). In practice, both operations are performed simul- 


taneously. 
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2.h Linearization of MSS Collinearity Equations 
The equations represented by (2.18) and (2.19) are nonlinear 
in nature. In order to investigate image and object point coordinate 
deformations, and to perform space resection utilizing the method 
of least squares, it is desirable to generate linear approximations 
for them. The usual technique applied in photogrammetry utilizes 
a Taylor’s series expansion about some approximations for the variables 
involved. The resulting linearized equations, as customarily written 
in adjustments, will be of the form 

F = _X_ + _B_ _A_ + _F^ = (2.20) 

2,1 2,2 2,1 2 ,q q,l 2,1 2,1 

in which A represents the Jacobian matrix of the functions with respect 
to observed quantities, and V the vector of observational residuals. 

The B matrix in Equation (2.20) is the Jacobian of the functions 
with respect to the parameters of the transformation, and represents 
corrections to approximations for these parameters. The F c ’ vector 
in Equation (2.20) represents the numerical values of the functions 
evaluated at initial observation values and parameter approximations. 
The final estimates of the observations are given by 

= + 

2,1 


L 

2,1 


V 

2,1 


( 2 . 21 ) 
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in which L° is a vector of a priori values of the observations. 
For the problem under consideration. 


L 

2,1 



( 2 . 22 ) 


in which the image space coordinates for each point will be considered 
as observations. The final parameter estimates are 


X - + _A_ (2.23) 

<1*1 <1*1 q*i 


in which X° is a vector of approximations and q is the total 
number of parameters. Of these q total parameters, 3 will represent 
the object space coordinates (X, Y, Z), and the remaining p - q - 3 
parameters will represent those parameters necessary to the functions 
used to model the six exterior orientation elements . 

Consider, first the linearization of the collinearity equations 
as given in Equations (2.l8). The functions may be written in the 
short form 

F = U = 0 

1 W 

Fp = tan y/c + V = 

W 


0 


(2.24) 
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where 


pi 

V 


X - X “ 

= M 

c 

Y - Y_ 


3,3 

c 

_w_ 


Z - Z 
L- C— 1 


(2.25) 


The elements of exterior orientation are functions of time (and 
hence x) , and may he written in functional form as 


x c = x c (P X > x > 

Y c - Y c (P y , x) 

Z C = z c , 

(2.26) 

m = 0 ) (P , x) 

U) 

♦ -*(£*. *) 

K = K (P^, X) 

where P-^, Py, P^, P^ , P^, P^ are vectors of the specific parameters 
in each of the given functions. If polynomial forms are assumed 
for these functions, these vectors, will contain the polynomial 
coefficients of the functions. For example, if it is assumed that 
the aircraft altitude Z c is represented hy the second order polynomial 

Z c = P 0Z + P 1Z x + P 2Z x2 
then the P_ vector would he of order three and of the form 


X Z 

371 


p 0Z 

P 1Z 

l_P2zJ 
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The linearized form of equations (2.2U) may he written in matrix 
notation as 


F = F° + A 1 fiy + B A + Bp A 


2,1 

in which 


A* 


2,1 2,1 




2,3 3,1 2T^^Tl 2,1 


( 2 . 27 ) 


3,1 


A*o 

57T 


[6X 6Y «Z] 


[6X C 6Y c 6Z c fiw 6<t> 6<] 


( 2 . 28 ) 

(2.29) 


and A-p B and B 0 are Jacobian matrices of the collinearity equations 
with respect to y, the (x, Y, Z) coordinates, and the exterior 
orient it ion elements respectively. The functional dependency of 
the Ap vector is now introduced, in linearized form; 


sx c 

= 

c x 

l,r 

J*x 
r ,1 

+ 

fix 

6* c 

= 

A 

l.B 

S,1 

+ 

< 3 y 6X 

6Z C 

= 

°z 

l,t 

*z 

t,l 

+ 

a z sx 


s 

C U> 

*(0 

+ 

Am «* 

fi<}> 

= 

A. 

l,g 

A 

g»i 

+ 

^4, 

fi* 

= 

C K 

l,h 

h,l 

+ 

d* sx 


(2.30) 
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in which the £ terms are Jacohian matrices of the functions assumed 
with respect to the parameters, the vectors are corrections to 
the functional parameters, and the d terms are derivatives with 
respect to the observed quantity x. Using (2.30), equation (2.27) 
becomes 




(2.31) 




Ill 


To get (2.32) in the conventional form of (2.20) then 



where p = r+ s + t+ f+ g + h = the total number of parameters 
in functions of the exterior orientation elements. 

In order to evaluate the matrices in Equation (2.32), the 
following relationship for a quotient of two functions, may be utilized 


H 3G - G 3H 


3 “g" 

= 

®P- __ 

3p - 1 Qg 

- G 3H 

(2.33) 

3p L_H 



H \jp 

H 3p__ 



Where p is any independent parameter. 
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Utilizing this relationship the Jacobian matrices may be evaluated: 


B 

o_ 

2,6 


w 

c 

9F 

9Y 1 

c 

9F 

9Z 1 

9F 

9uT 

9F 

9<l> 

9F 

9k 1 

3Fp 
9 X 

c 

Mp 

9Y 

c 

9Z 2 

c 

Mp 

9w 

^2 

3<j) 

Me 
9 k 


1 

W 


K + W m 31 3 [ " m 12 + W n 32 ] . 


^” m 21 + ¥ m 3 l ^ ^" m 22 + ¥ m 32 ^ 


^"“13 + W m 33 -* 


^” m 23 + W m 33 -* 


j [m 12 ( Z-Z c ) -m 13 ( Y-Y c ) -U/Wfa^ Z-Z c } “ m 33 ( Y ” Y c } } ] 


1 [m 22 ( Z “ Z C 5 " m 23 ( Y " Y c ) “ V / W{m 3 2 ( Z " Z c J ” m 33 ( Y “ Y c } } 3 


[-¥ cos k -U/¥ (U cos k -V sin tc) ] [V] 

[ W sin k -V/¥ (U cos k -V sin k)] [ — U ] 


-~L 

2,1 


~9F “ 


0 

9y 


2 



1_ sec (y/c) 

L^y J 


c 





(2.3*0 


(2.35) 
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B 

2,3 


¥ 

3F, 

3Y 


W 

1*2 

8Y 

1*2 

3Z 


(m i;L - U m 31 ) (m 12 - U U 13 “ g * 33 * 
xx w W W 

(mg! - I m 31 ) (m22 - I^) (*23 ~ 1 m 33 ^ 


(2.36) 


The C, D matrices will depend upon the functional forms assumed 
to model the elements of exterior orientation. As an example, if 
second order polynomials are assumed for these functions , then 


1 x x 2 0 0 0 . 

0 0 0 1 x x 2 


C = 

6 ,18 



0 0 0 


1 X x 2 


PlX + 2 ?2X x 
P1Y + 2 P2Y x 

D 

^1 


(2.37) 


(2.38) 


U?l*c + 2 P2 k x J 



llU 


Finally, the F° matrix may he written as 



U° 

W° 


tan (y°/c) + V° 
W°„ 


(2.39) 


in which the non indicates the evaluation at the original observa- 
tions and parameter approximations. 

An alternate linearized form may be obtained from Equation (2.19): 


F x = (X c - X) + (Z - Z ) U' - 0 

W 

F 2 = (y c - Y) + (Z - z c ) V’ - 0 

W' 


in which the auxiliaries U’ , V’ , W' are defined as 


~u r 


“ 0 “ 

V' 

= M t 

tan (y/c) 


3,3 


_ W L 


-3 


(2.40) 


For this case, the linearized equations may be written in the same 
form as Equation (2.32), and the Jacobian matrices may be evaluated 
as 
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_ 3F 1 

Ml 

IZa 

Mi 

M j. 

Mil 

3X C 

3Y c 

3Z 0 

3u) 

3* 

3k 

iFo 

L 3X c 

3£, 

3Y 2 

c 

Mg 

|2e 

3w 

Mg 

a<(> 

SF p 

3k_ 

-U' 


- U'V'(Z - 

z e ) 

9 

( 

1 

W 1 


w 

\Z 


1 

i 

-V’ 

_ 

“l + 

V' 2 “| 

( z - 

z e ) : 

W’ 


L 

W’ 2 J 

C i 

\ 
i 


Z - Z (- V’ sinw + W' cosu + U’ 2 costo) 
W' W' 

Z - Z rt (U' sinw + V* U* cosw) 

“IT “ 


tan(y/e)(-nim + U’ ,in n ^ ) Z - Z ( 

11 w> 3 ~w^ 

tan(y/c)(-nL + V' m ) Z - Z 
^ W’ 13 V 


(2.U1) 
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And the D matrices would be as previously defined. The JF T ° 
vector is 


F’° 

2A 


(X c ° - x°) + (z° 
(Y° c - Y°) + (Z° 


z° J 

U'° 

c' 

W’° 

c 

o 

IS3 

ytO 


W'° 


(2.U4) 


These linearized collinearity equations may also be used to derive 
the differential displacement equations (2.16). (see Appendix A). 

It has been the purpose of this chapter to provide the deriva- 
tion of the most general forms of analytic expressions which may 
be used to analyze the geometric aspects of MSS data recording. 

In the following chapter these general expressions will be used to gen- 
erate algorithms for the purpose of restituting MSS data arrays. 

In addition, alternate restitution methods based upon stochastic 
concepts rather than the known geometric relationships will be 
introduced. 



3. RESTITUTION OF MSS DATA ARRAYS FROM SINGLY SCANNED AREAS 


The term "restitution", as used in this investigation, indicates 
the processing of MSS data arrays in an attempt to generate within 
the spectral arrays, element positions which will represent a plan- 
imetric orthographic projection. If a display of the arrays is 
done after such processing, the resulting image would appear as 
a map substitute, much like the more conventional orthophoto. 

The term "rectification" has also been used by other investigators 
to denote this process. 

In this chapter, methods of restitution for ICS digital data 
arrays are presented, and the advantages and disadvantages of each 
method enumerated. In all of the discussion presented, it is assumed 
that only data from singly scanned areas are available. In the first 
section, procedures utilizing resampling algorithms are presented. 
These algorithms are useful as an initial step in order to simplify 
subsequent, more refined procedures. The second section deals 
with parametric methods using collinearity equations and polynomials 
which are derived from the linearized form of collinearity equations. 
Alternative methods based upon stochastic rather than deterministic 
concepts are presented in section three. These methods are denoted 
as nonparametric methods. Section four contains a discussion of 
the problem of assigning element elevations. These elevations may 
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be utilized in the parametric methods of restitution. Finally, - 
special techniques useful for restitution of MSS digital data, as 
opposed to continuously recorded data, are presented in section 
five. 


3.1 Resampling Algorithm 

Under the assumption of ideal orientation of the sensor, the 
data generated by the MSS in the direction of flight , or X direction , 
represents an orthographic projection with no further processing. 
Further processing would be required, however, in the direction 
perpendicular to flight, along each scan line, to remove the scan 
angle effects and topographic relief displacement. This may be ac- 
complished using resampling algorithms which operate on a single 
scan line at a time. The purpose of such algorithms are to digitally 
resample each data array line such that each element in the resampled 
line contains spectral values representing ground elements having 
equal dimensions in the Y direction on a datum. After applying 
such algorithms , the resulting arrays may still be in error due 
to the following: 

1. ) The use of such processing is based upon the assumption 

of perfect orientation, and the effects of changes with 
respect to time in the orientation elements X c , Y c , Z Q , 
a), 4> , k are neglected, which may result in appreciable 
image position errors. 

2, ) A lack of perfect synchronization between the aircraft 

velocity and angular scanning velocity, as well as the 
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difference between the angular resolution size in the 
x and y directions (6 and y, see Figure l.l6), results 
in a differential scale in the x and y directions after 
processing. 

Nevertheless , these resampling algorithms have several advantages : 

1. ) They are computationally efficient. The processing of one 

scan line at a time requires only a small image point 
buffer within the computer, and the computing time for 
such algorithms is relatively small., 

2. ) The algorithms result in the correction of image displacements 

due to panoramic recording. The arrays resulting after 
processing "look better” upon display, and may thus be 
useful for interpretation and identification of control 
points . 

3. ) These algorithms may be useful as an initial step in process- 

ing. Subsequent more refined algorithms for further geo- 
metric processing may therefore be simplified. For example , 
the piecewise polynomial to be discussed in Section 3.2.2 
will contain fewer terms if the arrays are first processed 
in this manner. Some of the nonparametric techniques to 
be discussed in Section 3.3 require that the arrays be 
first transformed to a trend surface, and the resampling 
algorithms accomplish this. 


This section is devoted to the derivation of resampling algorithms . 
The derivations consider two cases, the first of which deals with 
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the correction of scan angle effects only, under the assumption 
that the object space points have a constant elevation, and the 
second deals with a generalization to include topographic relief 
displacements . 

3.1.1 Resampling to Correct Seen Angle Effects [95] 

The purpose of this algorithm is to resample the data so that, 
after resampling, every sample element in the data array represents 
an element with equal dimension in the Y, or along scan, direction. 

That is , every scan line after resampling contains the same number 
of elements as before resampling, but object space elements are 
equally spaced along the datum. The algorithm developed by Phillips 
[95] at LARS is presented here. 

Figures 3.1a and 3.1b represent the geometry of this resampling 
procedure, with 3.1b being an enlarged view in the vicinity of the 
sensor. Referring to these figures, the following parameters may 
be defined: 

h - flying height above datum defined as before by h = (Z c - Z), 
where Z is the elevation (assumed constant) of the datum 
surface 

W - length of the scan line in object space 

N - total number of samples per scan, constrained to remain 

the same before and after resampling 
n - number of samples to right of nadir (if symmetric scanning 
is employed, as is usual, then n = N/2 and ny = a) 

J - sample number after resampling 
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S'j - spectral value to "be placed in the J-th position by- 
resampling 

Uj - sample number in the same line in the array before resampling 
corresponding to the J-th resampled position, to be 
solved for by the resampling algorithm 
Lj - integer portion of Uj 
Cj - fractional portion of 

- spectral value stored in position Lj in the array before 
resampling 

dj - distance along datum from beginning point of sampling 
to position J in the resampled array 
The variables 0, y, a, and sensor constant c are as described 
in Section 1.4. 

With the above variables, several relationships may be written. 

The scan line length is given by 

W = h tan(ny) + h tan(Ny - ny) (3.1) 

After resampling, every sample along the datum has the same dimension 
in the Y direction along the scan, defined by 

Y * W/N (3.2) 

From Figure 3.1a, the distance dj may be written, in terms of sample 
widths as 

d, = J(W/N) - (1/2) (W/N) = (2J - 1) W/N 

3 2 


(3.3) 
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This distance nay also be written in terms of sample angle (y) as 

dj = h tan(ny) + h tan(UjY - ny) (3.U) 

Equating the right hand sides of (3.3) and (3.M, and solving for 
Uj yields 

] (3.5) 


= ny + tan*" 1 

W 

hU 

[2J-1 

2 

- tan ny 

Y 


If Equation (3.1) is substituted into (3.5) » the final expression is 


Uj = n + y tan” 1 [(tan ny + tan(Uy - ny)) 2N - tan ny] (3.6) 


( 2 , 1 - 1 ) 


Figures 3.2 and 3.3 depict this relationship in graphical form. 
Figure 3.2 shows a plot of Uj versus j for a real data case, with 
Y - 6 mrad, and N = 222. Figure 3.3 represents the image displacement 


U - j versus j . The slight lack of symmetry in this example is 
J 

due to the fact that, for this particular flight, n ? N/2, 

If the quantities S'^, S'g, . . . S'j, . . . S’jj (see Figure 3.1a) 

represent the spectral values to be stored after resampling, a 

linear interpolation scheme may be used to assign the S' values. 

J 

If the quantity Uj computed above is considered to be composed of 
an integer part Lj and a fractional part Cj such that 


u j " l j + c j 


(3.7) 


12 U 



Figure 3.2. Resampled Array Column Positions 


Column Position Displacement (Ui-j) 
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Figure 3.3. Column Position displacement Due to Resampling 
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then the S'^ values may he interpolated from the spectral values 
originally stored in the array before resampling. The spectral 
value assigned to the j-th position would then he 


S' = S + C (S - S. ) (3.8) 

J L j 3 Lj+1 Lj 

Alternatively, if a nearest neighbor approximation is to he 
used, as advocated by Markarian, et . al. [65], the resulting function 
would he 


■ {> 
t b Lj+l 


Resampling in this manner results in the removal of the pan- 
oramic distortion within the data arrays. This processing differs 
from the usual reduction of panoramic image coordinates to a plane 
"equivalent frame photograph" position as presented by Case [ill]. 

The effective plane photograph position by Case would by y f = c tan 6 

J 

in which the panoramic image position is projected to a plane tangent 
to the cylindrical recording surface. The result of resampling 
is to project the panoramic image onto a plane at some effective 
principal distance c ’ . The case under consideration is illustrated 
in Figure 3.1b, for the symmetric case (n - N/2) . Since the total 
number of samples is constrained to remain the same , and each sample 
will be displayed with the same dimension before and after resampling, 
the distance Ob before resampling and O'b* after resampling must 
be the same. This is equivalent to imaging on a plane, indicated 
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by the dashed line in Figure 3.1b, lying at some equivalent principal 
distance from the perspective center and given by 

£ 

c f = k a (3.10) 

where k a is a function of the number of samples along the line 
and the angular resolution (y). 

From Figure 3.1b 


hence 


O’b 1 = Ob 


(c/k 0 ) tana = ca (3.11) 

for symmetric scanning of a on either side of the scanning axis. 

The above equation may be solved for k 

a 

tana 

k a - a (3.12) 

and subsequent image positions on the resampling plane are given by 

y f = c' tan(UjY - a) = c' tanBj (3.13) 

A simplified example may be illustrative. Consider the case 
of symmetric scanning, in which n ” N/2 and ny = a* For this case, 
the resampling algorithm may be written 

1 _ [ ~2.1 - 1 

Uj = n + y tan 1 [_ N tan ny - tan ny 
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or 


1 f2,1 - 1 

Uj = n + y tan -1 [_ N tan y - tan y 

If it is assumed for this simplified example that the total number 1,5 
of samples per line is eight (N = 8), and that the angular sampling 
interval (y) is 0.1 radian, then the total scanning angle on either 
side of the scanning axis (a) will be 0.4 radian. Referring to 
Figure 3.1b, to determine the spectral values to be placed in the 
third sample position in the resampled line (j - 3), 

JL _ R 2)(3) ~ (1) 

= 4 + 0.1 tan 1 L_ 8 tan 0.4 - tan 0. 

or 

U 3 = 2.26 


This represents the sample position along the line in the uncorrected 
array which would be used to assign spectral values to the j - 3 
sample position in the resampled array. If the linear interpolation 
method is used, then 


C 3 = 0.26 

and the spectral value assigned would be given by 

S' 3 * s 2 + 0.26(S 3 - s 2 ) 



If, on the other hand, the nearest neighbor method is applied, 
then the spectral value assigned to this resampled position for 
each channel would be given by 


S- 3 = S 2 


3.1.2 Resampling to Include Topographic Displacement [ 96 ] 

The algorithm given by Equation (3.6) is valid only for the 
case in which element elevations are assumed constant , i.e. « % 
for all j . A generalization of this algorithm was developed by 
Trinder [ 96 ] and corrects simultaneously for scan angle effects 
and image displacement due to terrain relief. The result to follow 
is a slight modification to that algorithm. 

Figure 3.U shows the geometry of this resampling. The illustra 
tion is similar to that of Figure 3.1a except that the actual topo- 
graphic surface has replaced the assumed flat plane representing 
the terrain. It is assumed that an elevation value, , is 
available for each element in the uncorrected array. With this 
assumption, the variables given in the preceding derivation are 
used except that a flying height above terrain may now be defined 


for each element as 
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The derivation proceeds in a manner similar to that given for 
the preceding case, and the resulting algorithm may he given by 
(See Appendix B) 


I . f JL pLL^-i 

Uj - n + y tan M h L |_ 2N h. 

+ 2N tan(ny) - h^ tan(ny)J 


• N tan(Ny - ny) 
1 


(3.1*0 


The spectral values may be assigned after resampling using 
either of the two methods cited in the preceding section, i.e., 
by linear interpolation of spectral values from the uncorrected 
arrays, or by a nearest neighbor approximation. 


3.2 Parametric Methods of Restitution 
In this section, methods of restitution based upon the MSS 
collinearity equations and piecewise polynomials are discussed. 

These methods are called "parametric methods", since an attempt 
is made in them to estimate scanning parameters . In both of these 
methods, an attempt is made to estimate the parameters of functions 
assumed to represent the behaviour of the scanner orientation elements. 

3.^.1 Restitution by MSS Collinearity Equations 
Probably the most logical mathematical technique to generate 
restituted arrays of singly scanned data is the direct use of the 
MSS collinearity equations as given in Equations (2.19). 'TTiese 
equations represent most closely the actual physical situation as 
it occurs during recording. After the assumption of appropriate 
functions to represent the behavior of the exterior orientation 
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elements, a solution may "be obtained using sufficient object space 
control points to apply least squares. For this purpose, the ex- 
panded linearized form of the collinearity equations, as presented 
in Equation (2.32), may be utilized. An iterative technique employing 
successive relinearizations about updated estimates from each 
preceding iteration is utilized. The final result is the determination 
of other po?Lnts at any desired scale, as a restituted array for 
each channel representing a single wavelength band. The method 
is subject to several limiting factors . 

1. ; Functional forms, such as polynomials or harmonics, repre- 

senting exterior orientation variations must be assumed. 

If the actual stochastic variation of these elements is 
not reasonably approximated by these functional forms , 
serious errors may result. To partially overcome this 
difficulty the flight line may be divided into sections, 
each having a separate set of coefficients. 

2. ) If the flight line is treated in sections, it is necessary 

to introduce constraining conditions at the section inter- 
faces. This insures obtaining the same results for each 
element along an interface from either side of the inter- 
face. (Thus the least squares solution must consider 
tne problem of parameter constraints . ) Tnis problem is 
addressed in more detail in Section 3.2.3. 

3. ) Some exterior orientation elements are highly correlated, 

and therefore recovery of all such elements may not be 
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possible using singly scanned data arrays. A previous 
example demonstrated this correlation for the <J> and X c 
terms. Under the conditions of flat terrain it is not 
possible, using singly scanned arrays, to differentiate 
between the effect of a 5<f> change and a 6X C change. This 
result was verified in real data tests using arrays obtained 
over terrain having relief differences of up to 300 feet (90 m) . 
Referring to Figure 1.21a as a first approximation, the 
change in X object space coordinate may be given as 

6X = (Z c - Z)<j> 

But the effect of this change may not be distinguished 
from the effect of a 6X c change. Thus, if functional forms 
are assumed for the X c , Z Q terms, the <j> term may not be 
recovered. Konecny [ 7 ] and Derenyi [22] note this problem 
and suggest that recording of selected exterior orientation 
parameters to sufficient accuracy represents a possible, 
although expensive, solution. Kratky [21] notes the problem 
as well, but suggests that acceptable restitution of the 
image may be accomplished even though absolute values of 
the elements will be in doubt. 

H.) If topographic relief displacement is to be compensated 

for during restitution, some method of assigning an elevation 
to each data array element must be devised. Since only 
singly scanned data arrays are being considered, the source 
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of such elevation information must he external to the data 
arrays themselves, such as existing topographic maps. 

The problem is further complicated in that these element 
elevations must be assigned prior to processing for resti- 
tution, i.e. , they must be assigned to the uncorrected 
arrays with acceptable accuracy. This problem is dealt 
with in Section 3.4. 

In summary, although the direct use of the MSS collinearity 
equations is accompanied by some problems which must be considered, 
these equations represent the best obtainable mathematical model 
for scanned imagery, and must be considered as the basic restitution 
method on which other approximations are based. In the next section, 
the use of polynomials derived from the linearized collinearity 
equations is investigated as a restitution method. 

3.2.2 Use of Piecewise Polynomials 
Since the recovery of absolute orientation parameters, based 
upon the MSS collinearity equations and availability of control, 
is not possible for each individual scan line, the use of such 
collinearity equations requires an approximation. It may be reason- 
ably argued therefore that a restitution based upon another math- 
ematical model could result in final element positions having posi- 
tional accuracies approaching those obtainable using collinearity. 

One such possibility is the use of piecewise polynomials. Using 
appropriate simplifying assumptions, these polynomials may be derived 
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directly from the more rigorous collinearity equations and thus 
represent a logical step in a search for an alternative mathematical 
model. 

Two cases may he derived. In the first case, it is assumed 
that the data are panoramically recorded. In the second case, the 
data are assumed to he resampled such that they are effectively 
recorded upon a planar surface (i.e. an attempt has been made at 
eliminating the panoramic effect) . 

3. 2. 2.1 Polynomials for Panoramic Recording . For this case. 
Equations (2.16) may he utilized directly. The object space coor- 
dinates are given by combining the ideal positions from Equations 
(l.ll) and (1.12), and the differential changes from Equations (2.16). 
If h = Z c - Z, and both Z c and Z are assumed constant, then h is 
a constant, and the transformation from the two dimensional image 
space into an object space which has been approximated as two di- 
mensional may be written 


X - X c + 6X C - h5<j» - h tan 6 

Y = Y c + <5Y C 4 h tan 6 4 tan 6 6Z C 4 h(l 4 tan 2 0)6w 


(3.15) 


Any number of polynomial forms may be defined, depending upon 

the functions selected to represent the variation of exterior 

orientation elements within a section. As a first example the 

assumption is made that the change in each orientation element 

(i.e., 6X , 6Y , etc.) varies as a linear function of x within 
c c 


each section, or 
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«c = *X + B X x 

5Y C = Ay + By x 

«e = A z + B z x 

6 “ “ k u * B u x 

a* - x 

6k = \ + B K X 


(3.16) 


These relationships may Tbe substituted directly into Equations 
(3.15) vith the results, after regrouping, of 

X = (X c + A x - hA^) + (B x “ hB^)x — hA^ *fcan0 ~ hB^x t cm 0 

Y = (Y c + Ay + hAj + (By + hBjx + (h + A z )tan0 

+ x tan0 + hA^ tan 2 0 + hB w x tan 2 0 

If the first two terms of the series expansion are taken for 
tan0(tan0 = 0 + ^-0 3 ), a* 1 ** using 0 = y/c (c is the effective 

panoramic recording radius), the resulting polynomials are 

X = A x + AjjX + A 3 y + (A 3 /3c 2 )y 3 + A^xy + (A^/3c 2 )xy 3 ^ 

Y = + B 2 x + B 3 y + (B 3 /3c 2 )y 3 + B^xy + (Bj + /3c 2 )xy 3 

+ B^y 2 + (2B 5 /3c 2 )y 4 + (B,_/9c 4 )y 6 + Bgxy 2 
+ (2B 6 /3c 2 )xy 4 + (B 6 /9c 4 )xy 6 
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in which 

A 1 = X c + A X ~ hA * B 1 = Y c + A Y + hA w 

Ag a B x - hB^ B 2 ■ By + hB w 

Ag = -hA^/c B 3 = (h + Ag) / c 

Aj| = -hB K /c B^ = B^/c 

B 5 = h^/c* 

Bg = hB^/c^ 

The last six terms of the second equation of (3.17) will he zero 
if the system is roll stabilized (i.e. w * 0). Note that the original 
12 parameters of Equations (3*16) have been consolidated to only 
10 parameters, since the parameters for the <$X C and 6 4 > terms always 
appear together, and therefore cannot be distinguished using singly 
scanned imagery of flat terrain. These parameters are A^ (i = 1, U) 
and Bj (j = 1 , 6 ). 

If the exterior orientation variations within a section are 
assumed as second order polynomials of the form 


5X C = A x + B^x + CyX 2 

<5Y^ — Ay t ByX + CyX 2 

6Z c s h + B Z X + C Z x2 

6a) = A^ + B^x + C w x 2 

5 <#> = + B^x + C^x 2 

^ + B K x + C K x 2 


( 3 . 18 ) 


61c 
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and these relationships are substituted directly into (3.15) » the 
resulting polynomials are 


X = A 1 + + AgX 2 + A^y + (A^/3c 2 )y 3 + A 5 xy 

+ (Ac/3c 2 )xy 3 + Agx 2 y + (Ag/3c 2 )x 2 y 3 

5 ^ (3.19) 

Y * B 1 + B 2 x + B 3 X 2 + B^y + (B^/3c 2 )y 3 + B^xy 

+ (B 5 / 3 c 2 )xy 3 + B^y + (Bg/3c 2 )x 2 y 3 + Bjy 2 

+ (2B T /3c 2 )y 4 + (B ? /9c 4 )y 6 + Bgxy 2 + (2Bg/3c 2 )xy 4 

+ (Bg/9c 4 )xy 6 + B^y 2 + (2B 9 /3c 2 )x 2 y 4 + (B 9 /9c 4 )x 2 y 6 

in which 


Ai = X c + A^ — hA^ B^ 

Ag - B x - hB^ B 2 

A3 = C x - hC^ B 3 

A^ * -hA^/c B^ 

A^ * -hB < /c Btj 

A 6 = -hC K /c B 6 

b 7 

b 8 

B 9 


+ % + H. 

b y + 

Cy + hC u 
(h + A^)/c 

i-! 

B z /c 

C z /c 
hA^/c 2 
hB w /c 2 
hC^/c 2 


The last nine terms of the second equation of (3*19) will be zero 
if the system is roll stabilized (w = 0), since these coefficients 
are dependent upon to only. The original 18 parameters in 
Equations (3.18) have been reduced to only 15 in (3.19) because 
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of the inability to separate parameters for <5X C and 6<f> using singly 
scanned data from flat terrain. The parameters are A^ (i =* 1, 6) 
and Bj (j = 1, 9). 

In both of these examples it has been assumed that all orien- 
tation elements were modeled using the same order polynomial. 

This need not necessarily be the case, since different order poly- 
nomials may be used. (e.g. X c - A x + B x x * C x t 2 + D x x 3 , 

Yg — Ay ByX ■fr’ CyX 2 , Z c — A ^ ^ Bf^x, ~ A^) * 

These polynomials represent a transformation from a 2 dimensional 
image space into a 2 dimensional object space. The object space 
is actually 3 dimensional, and the assumption has been made that 
h is constant in order to perform this transformation. This assumption 
may be relaxed if Z coordinates have been assigned to the data arrays 
prior to processing. In this case, h is no longer a fixed constant, 
but instead, h = Z c - Z may be used in Equations (3.15), where 
Z denotes the elevation assigned to each point, and Z c is the constant 
nominal flying altitude taken from the flight log. For linear func- 
tions, the polynomials when Z is included take the form 

X * A 1 + AgX + Agy + (A 3 /3c z )y 3 + A^xy + (A^/3c z )xy 3 
- Z[- Aij - Agx + (A 3 /Zc)y + (A 3 /3c z Zc)y 3 
+ (Ai,/Zc)xy + ( Aj,/3c z Zc)xy 3 ] 

4 4 (3.20) 

Y - B 1 + B 2 x + B 3 y + (B 3 /3c z )y 3 + B^xy + (B^/3c z )xy 3 

+ B 5 y z + (2B 5 /3c z )y 4 + (B 5 /9c 4 )y 6 + Bgxy z 

+ (2B 6 /3c z )xy 4 + (Bg/9c 4 )xy 6 - Z[(B 5 c 2 /Zc) + (B 5 /Zc)y z 

+ (2B 5 /3c z Zc)y 4 + (B 5 /9c 4 Zc)y b + (B 6 c z /Zc)x + (Bg/Zc)xy z 

+ (2Bg/3c z Zc)xy 4 + (Bg/9c 4 Zc)xy b + (l/c)y + (l/3c 3 )y 3 ] 



in which 


iko 


A 1 

“ x c + A X ** z c A <f> 

II 

H 

FQ 

y c + Ay + : 

a 2 

“ B X “ Z c B <f> 

b 2 = 

By + 

a 3 

= -Z c A k /c 

B 3 " 

(z c + Aij)/' 

Alt 

* - z cV c 

Bit * 

B z /c 

A? 

= 

b 5 = 

Z C A t»)^ C ^ 


ii 

td 

-©- 

B g = 

z cV c2 


The Btj, Bg coefficients will become zero if the system is roll 
stabilized (w = 0) and the corresponding terms in Equations (3.20) 
will drop from the second expression. Similar expressions, with 
a corresponding increase in the number of terms within the polynomials 
may be derived under the assumption of second order polynomial 
variations for exterior orientation elements within sections. 

Again polynomials of different order could be assumed for the dif- 
ferent orientation elements. 

Other investigators, notably Kratky [21], have considered the 
possibility of approximating the behaviour of orientation elements 
within a section by harmonic functions. In this case, functions 
may be written for exterior orientation elements of the form 


SX c 

ss 


sin 

c X x 

+ 

B X 

COS 

C x x 

«c 

— 

a y 

sin 

CyX 

+ 

By 

cos 

CyX 

sz c 

= 

h 

sin 

c z x 

+ 

B z 

cos 

C z x 


= 

A U> 

sin 


+ 

B(i) 

cos 

C w x 

6<t> 

- 

% 

sin 


+ 

B * 

cos 

V 

Sk 

ss 

Ak 

sin 

C K x 

+ 


cos 

C K X 


(3.21) 



These expressions may be substituted into Equations (3.15), which 
results in a complicated harmonic equation set. Kratky, in a 
subsequent communication with this investigator, indicated that 
use of these harmonics was troublesome. It vas found that, unless 
very close estimates of the parameters in these functions were known, 
it was difficult to obtain convergence for the parameters using 
least squares. 

3. 2. 2. 2 Polynomials for Equivalent Rectilinear Recording . The 
polynomials in Section 3. 2. 2.1 may be simplified if resampling is 
first done using the resampling algorithms presented in Section 3.1. 

The resulting data arrays will represent the data form as if the 
arrays had been recorded on a planar surface at an effective principal 
distance of c* , as presented in Section 3.1.1. The algorithm of 
Equation (3.6) is first considered, in which the j term would cor- 
respond to the y' parameter of Equation (3.13), and Uj would 
represent the corresponding position in the arrays before resampling 
and hence define the scanning angle 0. Since the y’ and 6 values 
are taken with respect to the scan line center, 

0 * (Uj - n)y (3.22) 

W 

y* * N (j - n - 1/2) (3.23) 

If the assumption is made that symmetric scanning is employed, in 
which n = N/2, and that resampling takes place upon a planar surface 
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at an effective principal distance of c' = c/ko as shown in Figure 
3.1b, the scan line length is given by 

W = 2c’ tan ny (3.24) 

The resampling algorithm for symmetric scanning is given by 




, IlSL 

= n + 1/y tan” 1 |_ N 


- 1 ) 


tan ny - tan ny 


] 


Substituting this expression for Uj into Equation (3.22) yields 


0 * tan 


-1 


J 


(2.1 - 1) 


N 


tan ny - tan ny 


] 


(3.25) 


Solving Equation (3.23) for 3 and substituting the result into (3.25), 
yields 



(3.26) 


Taking the tangent of both sides of (3.2 6 ) and substituting the 
expression for W from (3.24) results in 

Z± 

= c* 


tane 


(3.27) 


Ih3 


If (3.2T) is used in Equations (3.15) the resulting differential 
equations "become 

X = X„ + 6X. - h6$ - (h/c')y' 6 k 

( 3 . 28 ) 

Y * Y c + 6Y c + (h/c')y’ + (y’/c f )«Z c + h[l + (y'/c’) z ]5u> 

Assuming h is constant , in order to accomplish a transformation 
from two dimensional image space into a two dimensional object space, 
and that linear functions are used for the orientation element 
variations, as given in (3»l6), the resulting polynomials for this 
case are 

X - An + Ao x + A, y’ + A), xy’ 

1 ^ 3 (3.29) 

Y = B 1 + B 2 x + B 3 y’ + xy’ + B ? y ,z + Bg xy’ 2 

in which 

A 1 = X c + A X ' h A 4> ®i « Y c + A Y + hA w 


Aq = 

B X 

-"B ( 

b 2 = 

B Y + h B u 

A 3 = 

-h 

V°* 

B 3 " 

(h + A z ) /c 1 

A k = 

-h 

B k /c- 

B U = 

B z /c' 




b 5 = 

h kjc' 2 




Be * 

h B^/c 


If the scanner used is roll stabilized, the last two terms of the 
second equation of (3.29) will disappear. As mentioned before, the 
parameters associated with the 6<j> and SX C functions always appear 
together combined in the A^, A^ polynomial coefficients, and hence 
the effects of the 6X C and 6<|> terms are not separable, resulting 


Ikk 


in the reduction of the original 12 parameters to 10 for flat terrain. 

These parameters are A. (i = 1, k) and B. (j « 1, 6). 

•*“ J 

Note that the formulation of polynomials utilizing the resampling 
algorithm does not reduce the number of parameters . The parameters 
within these polynomials are dependent only upon the functions assumed 
for the elements of exterior orientation, and hence are unaffected 
"by resampling. The effect of resampling is to reduce the number 
of terms within the polynomials by making unnecessary an expansion 
to approximate the trigonometric function within Equation (3.15)* 
Assuming second order polynomials to represent orientation 
variations within sections as given in Equations (3.18) results in 

• X = Ai + iU x + Ao x 2 + Ai, y' + A, xy' + A/ x 2 y' 

J 5 (3.30) 

Y = + Bg x + x 2 + B^ y' + B,. xy' + Bg x 2 y f 

+ Bj y’ 2 + Bg xy' 2 + B^ x 2 y' 2 

in which 




X + A v - h A,* 
C x '9 

h. 

— 

Y C + A Y + h A u> 

h. 

= 

B x - h B * 

B 2 

= 

By +hB, 

A 3 

= 

c x - h °+ 

B 3 

S 

Cy + b C w 

A U 

s 

-h A^/c* 


- 

(h + A z )/c’ 

*5 

= 

~h B (C /c t 

B 5 

= 

B z /C ’ 

*6 

= 

-h C^/c* 

b 6 

= 

C z /c- 


B 7 = h Ajc' 

B e = h V c ’ 2 

B 9 = h Cjc<* 
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Similar to the presentation in Section 3. 2. 2.1, the effects of 
topography may he included within these polynomials by relaxing the 
assumption of constant h and including the element elevations assigned 
(h = Z Q ~ Z) for each element. Under the assumption that orientation 
element variations are linear, the resulting polynomials become 


X = A^ + A^ x + Ag y' + A^ xy' - Z[- Atj - Ag x 
+ (A 3 /Z c )y* + (A^/Z c )xy ' ] 

Y = B-l + B p x + B 3 y» + xy' + B ? y' 2 + Bg xy’ 2 
- Z[(B 5 c t2 /Z c ) + (B 5 /Z c )y' 2 + (Bg c ,2 /Z c )x 
+ (B 6 /Z c )xy» 2 + (l/c T )y f ] 

in which 


(3.31) 


A 1 = X c +A X~ Z c A * 
h. ~ B x - z c B <|> 
a 3 = “ z c Ac/® 1 

A^ = - Z Q B < /c f 


A 5 = A $ 


B 1 = Y c + A Y + Z c A w 
B 2 = B y + Z c B^ 

b 3 = < z c + V /c ’ 

\ - V c ’ 

b 5 - Z c V c ' 2 


H = B * 


B, 


" Z cV' T: 


and the B^, Bg coefficients will become zero if the system is roll 
stabilized (o> = 0) . 

If second order variations in orientation elements are assumed, 
the resulting polynomial functions become 


= A^ + Ag x + Ag x 2 + A^ y * + A«j xy ' + Ag xV 1 
- Z[- A^ - Ag x ~ A^ x 2 + (A lt /Z c )y* + (A^/Z^xy' 
+ ( Ag/ Z c ) x 2 y ' ] 


X 
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Y « Bj + Bg x + x 2 + y* + xy f + Bg x?y f + y ' 2 

+ Bg xy ' 2 + B 9 x 2 y ’ 2 - Z[(B c , 2 /Z c ) + (B 7 /Z c )y ’ 2 
+ (Bg c , 2 /Z c )x + (Bg/Z c )xy f2 + (B^ c f 2 /Z c )x 2 
+ (B 9 /Z c )x 2 y ' 2 + (1/c ’)y f ] (3.32) 

in which 


A 1 


x e + % - z c \ 

B i 


+ % + z c 

A 2 

s 

B X ~ Z c B * 

B2 

= 

b Y + z o B U 

A 3 

- 

c x - z 0 c t 

B 3 

= 

^ + Z c c a 

a 4 

s 

- z c V c ' 

B l» 

= 

(Z c + A z )/o' 

A 5 

= 

- z c V 0 ’ 

B 5 

= 

V c ' 

A6 

= 

- z c ° K / 0 ' 

b 6 

= 

Cz/C 

A 7 

= 


Bj 

= 

z o V 0 ’ 2 

a 8 

* 

B * 

b 8 

= 

z c V=' 2 

A9 

= 


B 9 

= 

z o V=' 2 

Bj , 

b 8 

, B 9 coefficients 

become zero 

if (D s 0. 


3.2.3 Constraints at Section Boundaries 
It has been stated previously that if a restitution is to be 
attempted using the collinearity or polynomial models , it is desirable 
to treat the flight line in sections. Dividing the flight line into 
sections allows the use of simpler functional assumptions for the 
exterior orientation elements to approximate their actual stochastic 
behavior with acceptable accuracy. If the entire flight line is 
treated as one unit, the functional behaviour assumed for the exterior 
orientation elements must be correspondingly more complex to yield 
comparable results. On the other hand, if the flight line is treated 


ll*7 

in sections , control points within each section may be used to determine 
functional parameters which are valid only within that section. 

This is desirable in an element by element restitution, since fewer 
terms must be evaluated within the computer for each element than 
if more complicated functional forms were assumed in treating the 
entire flight line as a unit. 

Treatment of the flight line in sections, however, introduces 
some problems. If each section is treated entirely independently, 
with an independent set of parameters, then restituted X, Y coordinate 
discontinuities will occur at the section boundaries . To avoid 
these discontinuities, introduction of constraints at section boundaries 
is necessary. 

In using collinearity equations or polynomials, the arrays 
are segmented along lines of constant x coordinate, i.e., constant 
data array line number. The joining conditions consist of the 
requirement that restituted positions along the line have the same 
values after transformation by either of the two possible collinearity 
or polynomial functions. Other constraints may also be considered 
in which, not only positions, but derivatives of first, second, 
or greater order may also be constrained to be equal along these 
section borders and thus a smoothing of the imagery across section 
boundaries may be introduced. These, however, are unduly complex 
and for this investigation, only the simplest case of constraints 
on coordinate position along the interface are considered. This 
problem may be efficiently treated if the concept of parameter con- 
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straints is applied at the time of estimating the parameters by 
least squares . In this case the linearized model may be -written 

AV + BA+BA. + F° = 0 

C A + G° = 0 (3.33) 

in which linearized constraint equations are used to augment the 
original linearized condition equations. 

If, for example, the polynomial forms given by Equations (3.20) 
are considered as a mathematical model, the constraint equations 
along the interface between one section, u, and the succeeding section, 
v = u + 1 may be determined as follows. 

Considering the first equation of (3.20) and equating X coordinates 
from sections u and v yields 

(A 1 + Ag x + y + (A 3 /3c 2 )y 3 + A^ xy + (A^/3c 2 )xy 3 

- Z[- A^ - Ag x + (A 3 /Z c )y + (A 3 /3c 2 Z c )y 3 
+ (Aj | /Z c )xy + (A^/3c 2 Z c )xy 3 ]} 

“ {A 1 + Ag x + A^ y + (A 3 /3c 2 )y 3 + A^ xy + (A it /3c 2 )xy 3 

- Z[- A 5 - A 6 x + (A 3 /Z c )y + (A 3 /3c 2 Z c |y 3 
+ (Ai + /Z (j )xy + (A] + /3c 2 Z c )xy 3 ]} y 
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Equating the constant terms and the coefficients of the y and 
Z terms respectively yields 

(A 1 + *2 x) u = (A 1 + *2 x) v 

(A 3 + x) u * (A 3 + A^ x) v (3.34) 

(A 5 + Ag x) u = (A 5 + Ag x)y 

Equating the coefficients of y 3 , yZ, and y 3 Z yields the following 
additional equations. 

[(A 3 /3c 2 ) + (Aj + /3c 2 )x] u = [(A 3 /3c 2 ) + (A^/3c 2 )x] v 

[A 3 /Z c ) + (A 4 /Z c )] u = E(A 3 /Z c ) + (A U /Z c )x] v (3.35) 

[A 3 /3c 2 Z c ) + (A^/3 c 2 Z c )x3 u = [(A 3 /3c 2 Z c ) + (A 4 /3c 2 Z c )x] y 

Since the c and Z c terms are treated as known constants for each 
flight line, it may "be seen that Equations 3.35 are redundant since 
each of them can he reduced to the second equation of (3.34). 

The remaining condition equations at each section boundary 
may be determined using the second equation of ( 3.20) , and enforcing 
the constraint condition of equal Y coordinate along the boundary. 

The condition equations which result by equating the constant terms, 
the coefficients of y, and the coefficients of y 2 in this constraint 
are as follows. 

(B x + B 2 x) u = (B^ + B 2 x) v 

(B 3 + B k x) u * (B 3 + B u x) v (3.36) 

(B 5 + Bg x) u = (B 5 + Bg x) v 
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The additional constraint equations which may be formed by equating 
y 3 , y 4 , y 6 » Z, Zy 2 , Zy 4 , Zy and Zy 3 coefficients, respectively, 
are 

[(B 3 /3c 2 ) + (B^/3c 2 )x3 u = [(B 3 /3c 2 ) + (B u /3c 2 )x] v 

[(2B 5 /3c 2 ) + (2B 6 /3c 2 )x] u = [(2B 5 /3c 2 ) + (2B 6 /3c 2 )x] v 

[(B 5 /9c 4 ) + (Bg/9c 4 )x] u = [B 5 /9c 4 ) + (Bg/9c 4 )x] v 

[(b 5 c 2 /z c ) + (b 6 c 2 /z c )x] u = [(b 5 c 2 /z c ) + (b 6 c 2 /z c )x3 v 

[(B 5 /Z c ) + (B^/Z c )x] u = [(B 5 /Z c ) + (B 6 /Z c )x] v (3.37) 

t(2B 5 /3c 2 Z c ) + (2B 6 /3c 2 Z c )x] u = [(2B 5 /3c 2 Z c ) + (2B 6 /3c 2 Z c )x] y 

{(B 5 /9e 4 Z c ) + (B 6 /9c 4 Z c )x] u = [(B 5 /9c 4 Z c ) + (B 6 /9c 4 Z c )x3 v 

1 /c * 1/c 

l/3e 3 ~ l/3c 3 

All of these equations are either identities, or may be reduced to 
one of the equations of (3.36) and are thus not independent. The 
final set of six independent equations necessary to enforce the 
constraint condition at section boundaries when utilizing the method 
of piecewise polynomials may be written by combining Equations (3.3*0 
and (3-36) as follows 

(A 1 + Ag x) u = (A 1 + Ag x) y 

(A 3 + A^ x) u = (A 3 + A^ x) v 

(A 5 + A gx) u = (A 5 + A g x) v (338) 

(Bi + B 2 x) u = (B^ + B 2 x) v 

(B 3 + B^ x) u = (B 3 + B^ x) v 

(B 5 + B 6 x) u = (B 5 + Bg x) y 
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The x image coordinate value in these equations represents 
the constant image x value at the houndary. Thus these equations 
contain only parameters and may "be linearized and used to augment 
the linearized condition equations during the process of least squares 
adjustment. Ix’ higher order polynomials are assumed for the elements 
of exterior orientation, the number of constraint equations will 
remain the same. However, each constraint will contain higher order 
terms in x. 

If the collinearity equations are utilized the necessary con- 
straints may be written quite simply by enforcing the conditions 
that , along the boundary between sections , the values of the exterior 
orientation elements must be the same whether computed using the 
parameters from section u cr those associated with section v. 


x cu 

= 

X C y 

Y cu 

= 

y cv 

z cu 

SB 

Z cv 

U) u 

= 


*u 

= 

*v 

K u 

= 

K v 
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For example, if linear variations are assumed for these elements 
within each section, the necessary constraint equations are 


( Ay + By x).^ 
(Ay + By x) u 
(Az + B Z x ^u 

(Ajj + B w x) u 

(A * + B * *) u 

(Ac + b k x)^ 


( Ay + x) y 
( Ay By x)y 

(Ag B^ x)y 
(\i * B u *>v 
+ *>v 

(Ac + B k x) v 


It is noted that the number of equations is the same as in the 
piecewise polynomial formulation, since the piecewise polynomials 
simply represent a different parameterization based upon a linearized 
form of the same collinearity equations. 


3.3 Nonparametric Techniques 

The preceding parametric methods of restitution attempt to 
estimate the coefficients of polynomials used to approximate the 
stochastic variation of scanner exterior orientation elements. 
Completely different techniques for restitution of remote sensing 
imagery have been applied by Leberl [45* 68], which may be termed 
"nonparametric” or "stochastic". The methods outlined by Leberl 
are in turn based upon earlier work by Pinkwart [97] and Schatz [ 98 ] - 
In this case, the difference in coordinate values in the x and y 
directions between uncorrected and restituted arrays is assumed to 
be a realization of a two dimensional stochastic field and attempts 
are made to obtain estimates for these differences by nonparametric 
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interpolate. ve methods . No effort is made vith these methods to 
estimate parameters associated with orientation. An estimate is 
made of the correlated portion of these coordinate differences, 
under the assumption that the mathematical model governing these 
correlated effects (often improperly called ’'systematic’' effects) 
is not known. 

The use of these methods assumes that a data array is a realization 
of a stationary random function, i.e., a random function for which 
the associated probability distribution functions remain unchanged 
with a shift along the independent variable axis or axes. Further, 
the property of ergodicity is assumed for the data array, which 
allows the use of values at different points in a field from only 
one of its realizations instead of many values at one point, which 
would require many realizations of the random field. References [ 109 ] 
and [ll^] may be consulted for a detailed discussion of stationary 
random functions and ergodicity. 

In addition, for certain of the nonparametric techniques to be 
discussed, particularly the linear ones, it is desirable to assume 
that the data array is isotropic. That is, that the correlated x 
and y mismatches between uncorrected and restituted data arrays 
are independent of coordinate direction, and dependent only upon 
distance between points. Thus, the field is assumed to be point 
symmetric. This assumption is definitely not valid for uncorrected 
MSS digital data arrays, and preprocessing is necessary in order 
to make this assumption applicable. 
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3.3*1 Arithmetic Mean [U 5 , 98] 

In this method, estimates of mismatch in the x and y directions 
are computed for each data array element based upon a known vector 
of such differences for a control point set. The control point 
vector of differences is defined as 


s xi ” x i “ x i 
s yi = Yi ~ Y i 


<3.1+0) 


for i * l,...n control points, where x^ , y^ represent the uncorrected 
image coordinate values for the i-th control point, taken from a 
gray scale display of the array, and are the corresponding 

object space coordinates of the same control point, after a conformal 
similarity transformation to the image coordinate system. Estimates 
of the x, y differences for an arbitrary element within the array 
are defined by a simple weighted arithmetic mean of the form 


l Vi s xi 
i-1 - 

n 

E Vi 
i=l 


n 

E Wi s yi 

i~l 

n 

E Wi 
i=l 


(3.U1) 


in which the Vi terms are weighting factors assumed to be dependent 
only upon the distance within the array between the element in 
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question and the control point being considered. Empirical weighting 
functions are then assumed [U5, 98] of tht? form 


or 



w. 


1 + dj® 


(3.U2) 


in which 

d^ = /(x - x^ ) ^ + (y- y i ) 2 


is the distance considered. The choice of the power m to be used 
is relatively arbitrary, depending upon the desired amount of cor- 
relation between widely separated points. If it is desired that 
correlation between points decrease rapidly, a choice of m = 3 
is suitable, in that correlation will drop off rapidly with distance, 
and only a nominal amount of computation is required. 

If this technique is applied directly using uncorrected data arrays 
and object space control which has been scaled to the image, the 
results will be acceptable in the x direction, but serious errors 
will remain in the y coordinate direction. The reasons for this are 
two- fold. 

1.) The arithmetic mean is a linear algorithm, both in x and 

y. While linearity of the data arrays can be well approximated 
over limited areas in the x direction, the recording of 
the data arrays in the y direction is panoramic in nature 


and hence non-linear. 
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2.) Due to the nature of scanning, the y direction in most 
instances will not have the same scale within the data 
arrays as in the x direction. (Scales in the x and y direc- 
tion will he the same only if the 0 and y sample angles 
are the same , and the period of rotation of the scanning 
prism is perfectly coordinated with the aircraft velocity.) 
Thus the data arrays , or images have differential scale 
values in the x and y directions. 

In order to effectively utilize the method of the arithmetic 
mean, then, it is desirable to pre-process the data arrays in a 
simple manner to circumvent the problems outlined above. To handle 
the first problem, each element y coordinate, or column number, 
would be processed using the resampling algorithm of Equation (3.13), 
in order to remove panoramic displacements. To resolve the second 


problem, these resulting column values would then be scaled to the 


object space control Y values, to remove differential scale. 


Thus , the processing of the data arrays by arithmetic mean 


would more realistically be defined by 


s xi = x i - x i 

s 'yi ” y 'i- Y i 


(3.U3) 


in which y'^ represents the image y coordinate values after pre' 


processing. 
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3.3.2 Moving Averages [45, 98] 

The method of moving averages represents a generalization 
which allows greater flexability in the point mismatch estimate 
computations than does the arithmetic mean. In using the method, 
the x, y mismatches are to be estimated using polynomials or other 
types of functions. The most commonly used examples are the linear 
(affine) case 


s x = a 0 + a x x + &2 y 
s y * b Q + b x x + b 2 y 


<3. MO 


or the second order polynomial case 


s x - a Q + a 1 x + ag y + a 3 xy + a^ x 2 + a^ y 2 
Sy “ b Q + ^ x + b 2 y + b 3 xy + b^ x 2 + b^ y 2 


(3.45) 


The a^, b^ coefficients of these functions are computed by 
weighted least squares utilizing the n reference points at which 
the x, y mismatches are known. The weight for each reference point 
is assigned as a function of the distance from the element where 
the estimate is being computed to the control (reference) point, 
such as the expressions of (3.42). For each point to be interpolated, 
a new set of the a^, b^ coefficients must be computed using 
weighted least squares. 

As an example, consider the affine case of Equations (3.44). 
Condition equations would be written for the control points where 
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x, y mismatch values are known. Considering the interpolation for 
the x coordinate, linear condition equations may be written of the 
form 



and the parameters (a 0 , a^, a^) may be determined by weighted least 
squares using the method of indirect observations (V + B i - F°). 
The weight matrix may be taken as 



An identical analysis and least squares solution technique may be 
used for the y coordinates to solve for the (b Q , b^ , bg) parameters. 
Having obtained a^, b^, Equation (3. WO is then used to compute 

A A 

s x , Sy. The entire procedure is repeated for every point to be 
interpolated. 

Note that as each x, y position is interpolated, the d^ terms 
of Equation (3.^2) will vary, hence the w^ terms will also vary. 
Therefore, a different coefficient set (a^, b^) must be evaluated 
for every point during a restitution, and the method may become 
computationally time consuming. 
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An advantage of the method lies in the fact that it will 
ordinarily require no pre-processing of the data arrays, as is nec- 
essary when using the arithmetic mean. The ’’systematic" effects 
of panoramic displacement and differential x, y scales are adequately 
incorporated within the estimating expressions. Such processing 
however, may allow the use of simpler functional forms to model the 
displacements in lieu of more complicated functions. 

It may be shown that the arithmetic mean method of the preceding 
Section (3.3.1) is a special case of the method of moving overages. 

In this instance, only the constant terms are used to estimate the 
mismatch at the point being interpolated, or 
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For this case, the linear condition equations become for the x 
coordinate 




l“ 


s xl 
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s xn 


from which, using the same diagonal weight matrix W x above, the 
parameter estimate becomes 



l6o 


n 

Z w i s xi 
i=l 

n 

Z w ± 
i=l 


Performing the same for y, leads to results identical to those of 
Equation (3.Ul). 

3.3.3 Meshvise L-inear Transformation [U 5 , 98 ] 

For this method, the control point set is first connected into 
a mesh. A commonly used mesh consists of a set of contiguous tri- 
angles. For each element to be transformed, this set of triangles 
is searched to determine the one in which the element being considered 
lies. Once this is determined, the transformation is accomplished 
by a linear affine transformation of the form given in (3.^U), 
in which the coefficients are determined uniquely from the three 
control points forming the triangle in which the element is found 
to lie. Three equation pairs of the form of ( 3 AM are written 
at the triangle vertices where coordinate displacements are known, 
and these are solved uniquely for the 6 coefficients (a^, b^) required. 

The method requires no pre-processing of the data, as the affine 
transformation allows for differential scale values in the x and y 
direction, and accomodates the panoramic effect by defining "secant 
planes" in the y direction. The method is computationally efficient 
in that only three points in the neighborhood of the element being 
transformed are used to uniquely define the coefficients, utilizing 
explicit formulas rather than least squares. The method gives 



identical results along the boundary between adjacent triangles 
when utilizing the parameters associated with either triangle. 

A variation of this method is termed overdeterained meshvise 
linear transformation. In this case, the control point field is 
used to define quadrangles, rather than a triangular mesh. For 
each quadrangle it is possible to define four possible triangles 
using three of the four possible points each time. Any element 
within a quadrangle will fall in two of these four triangles. 

The transformation outlined above may be performed on the element 
for each of these two triangles , and the final estimates of mismatch 
will be taken as the average of these two, or may be determined by 
the method of least squares. 

The method, however, suffers from serious drawbacks. The 
foremost of these is that it may not allow the transformation of 
all elements within the data array. In order that all elements 
within the data array may fall in one of the triangular meshes defined 
by the control points, it is necessary that, an a minimum, four 
of the control points fall at the extreme corners of the flight 
line. A control point, however, must be identifiable both on a 
map sheet (if maps are being used for planimetric object space 
control) and on a display of the data arrays. Clearly it is impos- 
sible that such points would fall at the flight line comers for 
every data array. Thus, the definable mesh using identifiable 
control points will not cover all of the flight line. This would 
make it necessary, if using a meshwise linear transformation, to 
define some element transformations by extrapolation utilizing the 
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parameters of Equation (3.M+), which in turn are based upon only 
3 control points. Such extrapolation is dangerous and may lead to 
gross errors, depending on the distance from the nearest control. 

A lesser, but still serious disadvantage of the method lies 
in the amount of human intervention necessary. The restitution 
of MSS digital data arrays, by whatever method, must be basically 
a computer oriented procedure to be compatible with automated inter- 
pretation when the two fields are finally combined. The necessity 
of defining triangular meshes by hand, which is necessary when using 
meshwise linear transformation, largely defeats the automated aspects 
so desirable in such arrays. 

3.3.^ Linear Least Squares Interpolation and Filtering 
Probably the most sophisticated of the nonparametric methods 
is that of linear least squares interpolation and filtering. 

The reader is referred to [99] for a detailed explanation of the 
method by Kraus and Mikhail. Basically, the method defines the x 
and y mismatch estimates in terms of observed mismatch values at 
the n control points and covariance matrices , written in the fora 

C " 1 

= 2,2n 2n,2n 2n,l (3.^7) 

2,1 

The vector of observed mismatches at control points, is considered 
to be composed of correlated random components , s_, and uncorrelated 
random components, r. 



im esaoxsn 
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1 = s + r 


(3.U8) 


The £ matrix is the autocovariance matrix for the n control points 
used, 


C 

2n,2n 


~xy 

n ,n n ,n 

5yx — yy 

n ,n n ,n _ 


( 3 . 1 + 9 ) 


and the £ matrix is the cross-covariance matrix of the element being 
transformed with respect to the control points. 


£ 

2n,2 





(3.50) 


The elements of both of these matrices are computed from 
covariance functions which are determined by a priori knowledge, 
or may be obtained from the control point data if a general form 
is assumed for the function, as described by Kraus and Mikhail [99]. 

If an isotropic field is assumed, then distances between points 
may be taken as the independent variable for generating the covariance 
functions. If this assumption is to be approximately valid, the 
data arrays must be pre-processed as described in Section 3.3.1, 
using the resampling algorithm of Equation (3.13) to reduce panoramic 
distortion, and scaled in the y direction. 


A more serious assumption associated with the method is station- 
arity. The method is predicated upon the fact that the random function 
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being considered for the field is stationary. If a non-stat ionary 
process is modeled as if it were stationary, then the computed sample 
covariance will not fit reasonably to an allowable covariance function 
form. One of the characteristics for an auto-covariance function 
associated with a stationary process is that its largest magnitude 
is at the origin. Thus, if sample auto-covariance values computed 
from the control point data increase in magnitude with increasing 
distance, this is evidence that the process is not stationary, and 
that one of the basic assumptions of the method of linear least 
squares interpolation is not being fulfilled. 

If this result is found within MSS data after pre-processing, 
it is possible to apply linear least squares only if stationarity 
is assumed for the stochastic field within a limited region about 
the point being interpolated. Therefore sample covariance points 
may be calculated only within a limited region, and allowable covar- 
iance functional forms fit to these points. The net result of this 
technique is that reference points beyond some critical radius about 
the point being interpolated are assumed to have insignificant cor- 
relation with the point , and thus do not contribute appreciably 
to the filtering and interpolation process. 

The theory of linear least squares interpolation and prediction 
has been used with success in geodesy, and its application in this 
field has been documented by Moritz [ 100 , 101 ]. In the latter paper, 
Moritz presents a generalization of the concept, in which correlation 
among points to be interpolated is considered. Leberl [1+5], however. 
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points out that work by Wolf [102] has shown that the inclusion 
of correlation between points to be interpolated leads to an over- 
constrained situation, and thus the independent interpolation of 
points as represented by Equation (3.1+7) represents the most feasible 
method of applying linear least squares. 

3. 4 Assignment of Element Elevations 
Element elevations must be pre-ass.igned to the uncorrected 
data arrays if the resampling algorithm of Section 3.1.2 is to be 
used, or if analysis or restitution is to be attempted using collin- 
earity (Section 3.2.1) or the more general polynomials (Section 3.2.2). 
The purpose of this section is to investigate methods by which this 
may be accomplished, and the problems associated with such methods. 

Element elevations may be derived from the MSS data itself 
only if the object space has been scanned more than once, as will 
be the case if overlapping strips of data are available. Most 
present data pertain to singly scanned areas, as were the data avail- 
able for this investigation. Therefore, elevations must be obtained 
from sources external to the MSS data. Possible sources for these 
data are existing topographic maps, such as the 1:21+000 scale topographic 
maps of the U. S. Geological Survey which were utilized for this 
investigation. For compatability with the MSS digital data, height in- 
formation must be processed into a digital terrain model (HTM). 

The DTM is a numeric representation of a topographic surface, in 
which elevations are given at discrete, closely spaced, X-Y point 
positions over the area being considered. 


The DTM may be of tvo types . 

1. ) Random, in which the planimetric point positions at vhich 

elevation values are known have a more or less random 
pattern over the area considered. This type of DTM would 
result, for example, if the digitization of a contour map 
was carried out by following selected contours, sampling 
X-Y point positions at intervals along each contour. 

2. ) Uniform pattern, such as a grid, in which the DTM is 

represented by closely spaced planimetric point positions 
arranged in a regular grid over the area. This type of 
DTM would result, for example, if the contour map were 
digitized by taking profiles at equally spaced intervals 
over the area. Then point elevations would be sampled 
at equally spaced intervals along each profile. 

In some cases, the DTM may be constructed from the random contour 
following method, and pre-processed into the patterned form. 

3A.1 Accuracy Restrictions 

At this point a method must be found to reliably assign an 
appropriate elevation to every data array element. A seeming paradox 
exists , since the DTM is essentially in an orthographic system, 
having been generated from map data. At the time when element ele- 
vations must be assigned, the MSS data arrays have all of the image 
displacements and distortions due to scan angle, topographic, and 
exterior orientation effects. However, the two data sets must be 
superimposed such that "acceptably accurate" element elevations 
may be obtained for the MSS data from the DTM. For purposes of 


this investigation, elevations were assumed to be "acceptably 
accurate" if the planimetric error due to height assignment would 
be less than a single resolution element. 

The solution to the problem of superimposing the two data sets 
lies in the nature of the ICS arrays themselves. Since the data 
arrays can resolve only to the value, ey, given in Equation (1.17) » 
it is only necessary to assign element elevations with an accuracy 
yielding a planimetric error of a magnitude equal to or less than 
ey, after subsequent processing. Recalling (1.17), the resolution 
in Y is 



y(Z c - Zj) sec 2 0 


Using Equation (1.19), the allowable height assignment tolerance 
Z e , is obtained as 



/tan 0 


y(Z. c - Zj) sec 2 0/tan 9 
y(Z c - Z^)/sin 9 cos 0 
2 y(Z c - Zj )/sin 20 


(3.51) 


Figure 3*5 shows graphically the magnitude of this allowable 
height assignment tolerance for various values of angular resolution, 

Y. Using Equation (3.51) at nadir, the allowable tolerance is infinite 
since a change in elevation has no effect on the image y position. 

The value of Z e decreases to a minimum at i*5° , where planimetric 
displacement becomes equal to height assignment tolerance. The 



Allowable Height Assignment Error (Z e ) as 
Percent of Flying Height Above Terrain 
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Figure 3.5. Allowable Height Assignment Error 
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allowable value then increases for scan angles beyond U 5 0 , since 
the ground size of a resolution element increases rapidly beyond 
this point. Such high values of scanning angle are rarely used 
in practice at the present time. 

To superimpose the two data sets within this tolerance, it is 
necessary to find a transformation which will yield a height assignment 
error of less than that given by Equation (3.5l). To do this, 
ground slopes must be considered. In Figure 3.6, the term SY represents 
the planimetric displacement in Y resulting from some transformation 
used to superimpose the MSS data upon the DTM for height assignment 
purposes. The e^ value is the size of the element at the given 
scan angle 0. The resulting elevation tolerance is given by 

Z e = X(6Y) (3.52) 

in which X is the ground slope in the vicinity. If the expression 
for allowable height assignment tolerance (ey/tan e) from Equation 
(3.51) is substituted for Z e into this expression and the resulting 
equation is solved for 6Y, then 

e Y 

SY = (3.53) 

X tan 0 

which represents the allowable planimetric tolerence which may be 
acceptable in an approximate transformation for height assignment 
purposes, and still result in less than a one resolution element 
planimetric error after subsequent processing, utilizing for example 
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the methods of Section 3.2. In reality, ground slopes rarely exceed 
0.5 and scan angles (e) normally do not exceed Uo°. Thus a trans- 
formation yielding errors on the order of 3 resolution elements 
may be appropriate for this step of height assignment . This value 
is necessary only near the flight line edges, and even larger toler- 
ances would be acceptable near the flight line center. 

A simple transformation which may be used for this step in 
the processing is the well known affine transformation. If this 
transformation is utilized within sections down the flight strips, 
the resulting planimetric error is found to be within acceptable 
limits. No constraints are necessary at section boundaries, since 
the transformation is used only internally within a program to 
approximately fit the data arrays to selected control points within 
the DTM. Coordinates of the control points required for this trans- 
formation may be taken from the same map source used to generate 
the DTM. The transformation is used only to assign element elevations , 
and no actual horizontal restitution of the M3S array positions 
would be done at this time. The more sophisticated approaches given 
in Section 3.2 would be utilized for the actual transformation of 
elements into a restituted array. 

3.H.2 Contour-Scan Intersections 

This method of element elevation assignment seems to be the 
simplest conceptual procedure if the DTM available is of the random 
type generated by the contour following procedure. The DTM data, 
if stored according to contours, may be used directly in the form 
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generated by the digitizing equipment. For this type of DTM, as 
mentioned previously, contours may be sampled in a digitizer at 
random intervals along their length. 

Each scan line may be transformed according to the affine 
transformation described in Section 3.^.1. The parameters of this 
affine transformation are estimated using the method of least squares , 
based upon ground control points. At this stage, the contours and 
the transformed scan lines will be superimposed approximately in 
a near-orthographic system. For each scan line thus transformed, 
all contours in the vicinity are searched digitally until an 
intersection of the contour and the scan line is detected. At this 
point, that position along the scan line is assigned the elevation 
of the contour. The tracing of the contour is then continued to 
determine if any further intersections occur with this same scan 
line. All contours in the area are traced in this manner until 
all intersections are found for this scan line. The intersection 
points are then re-sorted so that they are in order along the scan, 
and elements which are intermediate between intersection points 
are assigned elevations by linear interpolation between these known 
elevations. This process is repeated until every scan in the flight 
line has been processed. The resulting element elevations are stored 
on the original MSS data tape as an extra channel of data, in the 
same format as that given in Section 1.2.3. 

The advantages of the method are : 

1.) Relatively efficient computation: Only straight line 

equations are used throughout for finding contour-scan 
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intersection points and for interpolating element elevations 
between these intersections. 

2.) Small storage requirements: At any time, the computer 

is required to store only the positional data for a single 
scan line and a single contour. 

The major disadvantage of the method is the determination of 
the element elevations near the beginning and end of each scan. 
Contour-scan intersections will obviously occur exactly at the scan 
line ends only by chance, and with virtually zero probability. 

Thus , some method must be devised for assigning those elevations 
near scan line ends before the first intersection point and beyond 
the last. Extrapolation may be dangerous. An attempt must be made 
to interpolate these end points between contours, and a complex 
decision procedure is necessary to determine the contour pair between 
which an end point lies. 

3.^.3 Polyconic and Polynomial Surface Approximations 

A different approach to the problem of element elevation 
assignment is to attempt to adequately model by analytic expressions 
the topographic surface in the region of each element , based upon 
some subset of the digitized sample points . One example of a mathe- 
matical modeling technique is the multiquadric analysis advocated 
by Hardy [103, 10U]. The latter reference contains a n excellent 
summary and introduction to the method. 
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The general form of the mathematical model which Hardy advocates 
is 


n 

£ c k [q.(x k , y k , x, y)3 ■ Z (3.5*0 

k=l 

in which x k , y k > (k=l,...n) are the coordinates of the sample points be 
being considered in the vicinity of the element position x, y after 
its approximate transformation by the affine transformation, and 
z is the estimated element elevation associated with the point. 

The C k terms are coefficients to be multiplied by each of the quadric 
surface equations assumed. The C^, coefficients are estimated only 
once , and retain the same value over the section being considered 
regardless of the (x, y) location of the point being interpolated. 

The z elevation is a function of x and y resulting from a summation 
of a series of quadric surfaces having axes of symmetry at the x k , 
y k sample points. According to Hardy, the most often utilized forms 
of the quadric surfaces chosen are those of a cone, in which 

q(x k , y k , x, y) = [(x k - x) 2 + (y k - y) 2 ] 1 / 2 (3.55) 

of a hyperboloid, in which, 

q(x k , y k , x, y) = C(x k - x) 2 + (y k - y) 2 + C] 1 / 2 (3.56) 

of a paraboloid, in which, 

= [(x k - xj} 2 + (y k - y) z ] a ; a - 1/2 


<l(x k , y k , x, y) 


(3.57) 



175 


and an inverse hyperboloid, in which 

y k» x ’ = ^ x k " x ^ 2 + ^ y k “ y) 2 + C]~ l/2 (3.58) 

Thus , when given a set of n sample points such as would be 
generated, for example, in contour digitization, at every sample 
point an equation could be written of the form 

n 

Z C k [q(x k , y k , x ffi , y m )] = % (3.59) 

k=l 

in which m is the sample point being considered. There would be 
n such equations generated, assuming it is desired to have a quadric 
axis at each sample point. Using the matrix notation of Hardy 

A C = _Z_ (3.60) 

n,n n,l n,l 

the coefficient vector desired, C., could be solved for by 

C = A" 1 Z (3.61) 

n,l n,n n,l 

and elevations for elements could then be calculated from 
Equation (3.5*0. 

The system has many advantages, one of which is that it allows 
nonlinear modeling of the topography, which is obviously not possible 
with the linear interpolation scheme of the contour-scan intersection 
method. It also allows the computation of element elevations near 
scan line ends with more reliability than with contour-scan 


intersections . 
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However, if a large data set of points is available in the 
DTM, then the matrix A in Equations (3.60) and (3*6l) will become 
very large. The computing time necessary to solve (3.60) for the 
coefficients will be quite large, and the time necessary to compute 
each element elevation will also become excessive, as the summation 
of n polynomials to be evaluated must be performed. In cases of this 
nature, Leberl [**53 suggests using a relatively small subset of the 
total sample points available, such that 

d k * d 

where d defines some convenient "circle of influence" for each element , 
and d k is the distance from the element being considered to each 
sample point. 

^ = /(x k - x) z + (y k - y ) 2 

Another possible compromise to cut computing time would be to use 
a larger subset of the sample points (and hence larger A and C. 
matrices), but evaluate the multiquadric function of Equation (3.5*0 
only at selected element locations along each scan line, for example, 
every fifth element. Then points intermediate between these could 
be assigned elevations by linear interpolation. 

A recent paper by Jancaitis and Junkins [105] presents a 
polynomial procedure for modeling the topographic surface, with new 
polynomial coefficients being computed as the vicinity being considered 
is varied. In this method, the area is divided into a number of 
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small square areas to form a lattice. Within each square, the topo- 
graphic surface is to he represented by a polynomial of the form 

z = E Ch x 1 yJ (3.62) 

ij 

in which the coefficient subscripts correspond uniquely to the ex- 
ponents i, J of the x, y arguments. 

The coefficients of this polynomial, which ordinarily is 
chosen to contain 12 terms, are determined by matching the estimated 
elevations and slopes in the x and y directions of the four corner 
points for each square incremental area, thus enforcing continuity 
in elevation and slope at the boundaries between adjacent squares. 

These estimated elevations and slopes at each corner or lattice 
point are obtained by fitting a plane at each lattice point to some 
subset of the sample points determined during digitization, in the 
immediate vicinity of the lattice point. 

After the preliminary affine transformation of each element, 
these transformed coordinates are utilized directly in Equation 
(3.62) to obtain an estimate of the elevation for that element. 

3 .^.U Other Interpolative Methods 
The methods of Section 3.3 may also be utilized to estimate 
elevations of elements within the MBS data arrays. In this case, 
the interpolation becomes a one-dimensional case , in which the 
reference point set would be represented by the digitized values 


of the DTM. 
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The methods which may be applied may be taken directly from 
Section 3.3. For the arithmetic mean, for example, the estimate 
for elevation of an element would be given by 


n 


I v k z k 
k=l 


n 

S v k 


k-1 


in which the z k values are the elevations at the sample points of 
the DTM, z is the estimated elevation of the point being considered, 
and, the weighting factor w^ for each sample point may be evaluated 
by either 

w k = 1 or w k = 1 

d\ 1 + 

in which 

d k = /(x - x k ) 2 + (y - y k ) 2 

where x, y is the approximate orthographic element position after 
the affine transformation, and x k , y k is the digitized position of 
the point at which z k has been sampled. 

The DTM will generally contain a very great number of sample 
points. To obtain reasonable computing times it will normally 
be required that only points within some circle of influence will 
be considered, such that 

d k * d 


The same procedure as that outlined in Section 3.3.2 may he 
used for a moving average interpolation, vhich vill allow a non- 
linear interpolation. The estimate will, of course, be only a one 
dimensional estimate in z. The use of the meshwise linear transfor- 
mation is hardly practical for estimation of element elevations, 
in that it requires the determination of a unique set of contiguous 
triangular meshes in order to implement the method. With the very 
large data set of points generated in digitizing the DTM, the number 
of triangles becomes excessive. 

The use of linear least squares interpolation to estimate element 
elevations requires first the reduction to a reference, or trend 
surface. Leberl [1*5] points out that meaningful correlations between 
terrain points exist even over sizeable distances (30 km). He 
recommends the use of the method of moving averages to first estimate 
the trend surface, and subsequent utilization of the linear least 
squares, in an attempt to refine the element elevation estimate. 

Since, as was shown in Section 3.1*.l, the final positional accuracy is 
relatively insensitive to small errors in element elevations, it 
is doubtful that the increase in element elevation accuracy is worth 
the increased computing time necessary to implement the method 
of linear least squares interpolation. 

All of the above discussion has been concerned with general 
mathematical models which may be used for analysis or restitution 
of MSS data. In the following section, the unique problems associated 
with restituting MSS data stored in digital form are discussed. 
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3*5 Direct and Inverse Restitution Techniques 

The restitution equations described in Sections 3.2 and 3.3 
were based upon the premise that a transformation was being used 
from the uncorrected data array (or image) into the orthographic 
object space. These equation forms and relationships are most con- 
venient if estimation of the parameters of imaging is being attempted, 
and may be used with efficiency for restituting conventional recording 
media, such as photographic imaging. The use of such relationships 
may be termed a "direct" approach to restitution, since projective 
relationships are used to project directly from uncorrected x, y 
image positions to restituted X, Y positions (by analogue or analytical 
means ) . 

Because of the nature of digital MSS data, however, further 
considerations are necessary in attempting to apply this approach 
to such a digital data form. The digital data is stored in a line 
by line fashion, and within each line, spectral value locations are 
limited to integer positions, both in the uncorrected image array, 
and the restituted output array resulting from such processing. 

If the direct technique were used, then restitution would be 
accomplished by measuring x, y uncorrected image coordinates and 
projecting to restituted X, Y positions using, for example, the 
collinearity equations of (2.19) or one of the piecewise polynomial 
formulations. This leads to two problems. 

1.) The resulting X, Y position pairs will be randomly generated. 

This means that a large number of the coordinate pairs 
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must be stored in an output buffer. A continuous sorting 
must be implemented to arrange these pairs into positions 
along array lines to be stored in a restituted output 
array. 

2.) Unique X, Y pairs may not be generated for output. The 
X, Y pairs thus generated must be in integer form in order 
to occupy locations within the output array. The same 
integer position within the output array may result from 
two or more x, y positions in the uncorrected input array, 
and thus result in a multiplicity of spectral combinations 
which could be assigned to it. Further, there will be 
positions within the output array which would have no 
spectral values assigned because for some X, Y output array 
positions, there may be no x, y pair which would map into 
these positions through the projective expressions, since 
the X, Y values are limited to integer values. Spectral 
values for such "holes" may be assigned by interpolating’ 
spectral values from surrounding output buffer positions. 

The advantage of this method is that use of the direct projective 
restituting relationships (collinearity or polynomial) requires 
no iterative procedure. The element elevation values assigned are 
associated with the x, y uncorrected image positions. Hence, for 
each element to be projected, x, y and object space Z coordinates 
are known, along with the parameters of projection. The X, Y 
restituted position may be calculated with a single application 
of the restituting relationships. 
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An alternate restitution technique is to "step” along positions 
in the restituted output array as the independent variables (X, Y) 
as advocated by Markarian, et. al. [69], and compute the associated 
x, y position for each restituted array position. The spectral 
values in this position are then assigned to the restituted array 
position. For purposes of this discussion, such an approach will 
be called the "inverse" restitution technique. 

This procedure circumvents the storage and sorting of an output 
buffer within the computer, and assures that each point within the 
output array will have assigned to it a spectral value for each 
channel. 

However, this approach also has two major drawbacks associated 
with the data processing. 

1.) Since X, Y values are assigned as the independent variables, 
then x, y uncorrected array positions are randomly generated, 
thus necessitating the storage of a large number of these 
values in a buffer during restitution. After the computation 
of the x, y position for each associated X, Y output location, 
this buffer must be searched for the nearest integer array 
position to the desired x, y image position, and the associated 
spectral values stored in the restituted output array 
position. This input buffer, however, is somewhat easier 
to handle than an output buffer, since it is in array form, 
containing integer array positions in a line by line arrange- 


ment 
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2.) Perhaps more serious is the requirement of an iterative 
computational procedure. This may not apply for those 
based upon the greatest simplifying assumptions, such as 
the assumption of flat terrain and constant or linearly 
varying exterior orientation elements. The following 
analysis explains the procedure. 

For this purpose, the inverse form of the MSS collinearity 
equations (2.18) are considered. Capitalizing on the fact that the 
left hand side of the first equation of ( 2 . 18 ) is zero, these 
equations may be written as (dropping the j subscripts) 

0 * m 11 (X - X c ) + m 12 (Y - Y c ) + m 13 (Z - Z c ) 

(3.63) 

tan(y/c) = - [m^X - X c ) + Y - Y c ) + m^Z - Z c )] 
m 31^ x “* x c^ + m 32 (Y - Y C J + m 33 (Z - Z c ) 

Note that in the first equation of (3.63), the only variables are 
X, Y, Z, and x, which is included in the functions assumed for ex- 
terior orientation elements. Assuming for the present that Z has 
a constant value (perfectly flat terrain) then if X, Y are treated 
as independent variables, and if the parameters associated with 
the X c , Y c , Z c positional orientation elements and the w, <{>, k 
angular orientation elements are known from a previous space resection, 
then the first equation of (3*63) may be solved for x. This solution 
may be direct if fairly simple functional forms have been assumed 
for orientation elements, or iterative if more complex functional 
forms have been assumed. The value obtained for x may be substituted 
directly into the second equation of (3.63) to then solve for the 
y image coordinate. 
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If for this case the assumption of constant Z is relaxed, an 
iterative technique is required. If element elevations have been 
assigned to the data arrays before processing is begun, then they 
are associated with the x, y uncorrected image positions not the 
X, Y coordinate values used as the independent variables. An iterative 
technique would then be employed in which a Z value would be estimated 
(perhaps as an average elevation) and equations (3.63) would be 
solved for new x, y values. The element elevation stored in this 
location would then be taken as an updated value of Z, and new 
x, y values computed. This process would be repeated until no 
meaningful difference in the x, y values is noted. This is the 
technique to be utilized if a restitution based directly upon col- 
linearity is to be employed. 

If a polynomial approximation is to be used for an element 
by element transformation using this inverse technique, the problem 
becomes more complex. It is, in general, not possible to write 
a linearized transformation equation which may be solved explicitly 
for x, y image coordinates if X, Y object space coordinates are 
designated as independent variables, since the orientation elements 
of scanning are given as functions of the x image coordinates. 

If an element by element restitution is to be attempted by the 
inverse technique utilizing polynomial approximations, it is 
suggested that the polynomials of Section 3.2.2 be utilized in 
an iterative manner. 
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This choice between the direct and inverse methods of final 
transformation of elements is necessary only with the direct col- 
linearity and polynomial formulations. The resampling algorithms 
of Section 3.1 are already in the inverse form, in which the resti- 
tuted array position, j, is used as the independent variable, and 
the corresponding position, , in the original uncorrected array 
is computed in order to assign spectral values within the restituted 
array . 

The inverse technique is easily incorporated within the non- 
parametric methods. These techniques are not based upon analytic 
forms of modeling the transformation, but are based upon stochastic 
concepts in which an attempt is made to estimate the correlated 
components of a stochastic field without defining an analytic form. 
The mismatch estimates resulting from these methods may be applied 
to the X, Y positions to estimate the x, y uncorrected image position 
in order to assign spectral values. 

In summary, if direct restitution is attempted using projective 
collinearity or polynomial forms, one is confronted with problems 
of output buffer storage and sorting in generating restituted arrays, 
and the problem of assigning spectral values to the X, Y restituted 
position for which there are no associated x, y image positions 
based upon the direct projective relationships. On the other hand, 
if inverse relationships are used, the storage and sorting procedures 
are simplified, but the actual transformation calculations require 
an iterative procedure. For this investigation, no attempt was 
made to assess the relative efficiency of the direct approach with 
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respect to the inverse technique. The restitutions performed in 
this study utilized the inverse approach, in order that modifications 
of software routines already available at LARS could be used. It 
is left to future investigation to examine the direct technique. 

3.6 Concluding Remarks 

This chapter has served to introduce possible methods of 
restitution for singly scanned MSS digital data arrays , and some 
of the problems associated with such restitution methods. The 
resampling algorithms first presented are useful in many cases as 
a first step in the restitution process. The pre-processing of the 
data using these algorithms is useful in simplifying subsequent 
parametric algorithms, or in approximating isotropy for some of the 
nonparametric methods. 

The parametric methods of restitution utilizing MSS collinearity 
equations or piecewise polynomials represent perhaps the most 
intuitively obvious methods of restitution, since they are based 
directly upon the known relationships of the MSS transformation. 

If these methods are to be used to obtain the best possible accuracy, 
then Z coordinate information must be introduced into the data 
arrays from some source. 

The nonparametric methods included represent an entirely dif- 
ferent approach to the problem of restitution. These algorithms 
assume that no analytical form is known for a mathematical model. 

An attempt is made to estimate the correlated effects within the 
data arrays, based upon stochastic concepts. 
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The unique nature of the digital data form requires that special 
data handling techniques he utilized for an efficient element by 
element restitution. A choice is necessary between a direct approach 
utilizing the direct projective relationships relating the image 
and object space element locations, and an inverse approach based 
upon the imaging form of these relationships. The direct approach 
has the drawbacks of requiring the storage of a random output buffer 
coupled with sorting and interpolation requirements, but has the 
advantage that the direct projective relationships may be used and 
require only a single application for each element being transformed. 
The inverse approach results in simpler storage and searching proce- 
dures, but requires an iterative procedure for the element by element 
calculations required in the restitution. 

In the following chapter an operational procedure is presented 
to collect many of the concepts discussed previously into a logical 
sequence for analysis and restitution of digital data arrays. This 
procedure is used to test the methods of restitution on existing 
MSS digital data arrays, and the accuracies obtained from these methods 


are compared. 
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4 . NUMERICAL RESULTS 
4.1 The Data Arrays 

The methods of analysis and restitution presented in the pre- 
ceding chapter were tested on two flight line arrays of actual data 
available at LARS. Both were flown during the corn blight watch 
experiment administered by LARS during the summer of 1971 [106, 1073* 
The University of Michigan scanner was used at a nominal altitude 
of 5000 ft. (1524 m) above sea level, and digitization performed with 
a nominal angular resolution of 6 mi Hi radians. 

The first of these flights, number 208, was located in north- 
western Indiana near the city of Lafayette. Figure 4.1 depicts 
a grey scale display of this array from the line printer. The display 
is of channel 6, which represents the spectral band between the 
wavelengths of 0.58 and 0.65 micrometers. Figure 4.2 portrays this 
same flight line, also from channel 6, as an image taken from the 
digital display monitor. These images are at an approximate scale 
of 1:58000. This array contains 1591 lines and 222 columns, and 
represents a ground area approximately eight miles (13 km) long 
and one mile (1.6 km) wide. Since it was desired to utilize existing 
U. S. Geological Survey maps as a source of specific point positions, 
as well as terrain elevations, it was necessary to locate array 
positions on the line printer display which were also identifiable 
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on the map sheet. This flight line was well suited for this purpose. 
Terrain variation within the flight line was quite small, on the order 
of 50 ft. (15 m), ranging from an elevation of approximately 700 ft. 

(213 m) to approximately 750 ft. (229 m) above mean sea level. 

A large proportion of the area within the flight line was under 
cultivation, with no forest cover. The area contained numerous 
fence, road, and stream intersections which could be easily located 
on both the line printer display and the map sheet. These array 
positions were determinable on the line printer display to within 
one or two resolution elements. After extracting these array positions, 
the corresponding points were marked on the map sheet for later 
coordinate digitization. 

All points thus located on both the digital array and the map 
sheet are referred to in this section as reference points. From 
this set of reference points, a subset of control points was chosen. 
These control points served as the basic data source from which 
parameter estimates or nonparametric algorithms were calculated 
in the analysis methods to be discussed. The points in the remaining 
subset of reference points are referred to as check points, and 
are used to assess the accuracy obtained in the subsequent analysis 
methods. A total of one hundred such reference points were located 
within the flight line. Of these, thirty nine were used as the 
actual control points to determine the parameters for the various 
analysis methods, and sixty one were withheld as check points to 
investigate the results which were obtained. A mistake in digitizing 
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one of the check points was subsequently discovered, and the tests 
performed used only sixty check points. 

The second flight line chosen, flight number 218, is located 
in west central Indiana near the city of Bloomington. Figure 4*3 
shows a line printer display of this flight line, and Figure 4.4 
depicts the same flight line from the digital display monitor. 

As before, the images shown are for channel 6 representing the 
spectral wavelength band from 0.58 to 0.65 micrometers. These images 
are at an approximate scale of 1:54000. This array contains 1439 
lines and 222 columns, and represents a ground area approximately 
seven miles (11 km) in length, with a width of approximately one 
mile (1.6 km). Total terrain elevation variation within the flight 
line is on the order of 300 ft. (91 m) > ranging from a low elevation 
of approximately 550 ft. (167 m) to a high of about 850 ft. (259 m) 
above mean sea level. Most of the area is wooded such that a tree 
canopy obscures much of the ground. The area contains fewer roads 
and fences than the flight line previously discussed, with far 
fewer recognizable intersections and easily identifiable array 
positions. For these reasons greater difficulty was encountered 
in reliably locating reference points within the array. Thirty 
three points were finally assigned over the entire flight line, 
of which twenty four were utilized as control to determine ths 
parameters of the analysis, and nine were withheld as test points. 
However, because of discovery of a mistake in digitization of one 
of the control points, a set of only 32 was used (23 control, 9 
check) in performing the tests of analysis methods discussed below. 
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Uncertainties of 3 or 4 resolution elements were not unusual in 
assigning array positions to these 32 points. 

In order to illustrate the geometric distortions inherent in 
data arrays obtained by multispectral scanners, planimetry was traced 
from the map sheets for both flight lines for comparison purposes. 

Figure 4.5 shows some of the planimetric features of flight 208 
near the top of the flight line, and Figure 4.6 portrays planimetry 
from flight 218, near the bottom of the flight line. In these 
figures, a heavy solid line is used to denote roads, a lighter solid 
line is used to represent streams or ditches , and a light broken 
line represents fences. A comparison of these figures with those 
of Figures 4.2 and 4.4 will illustrate the deformations in the imagery 
which are perceivable to the human eye. The most readily apparent 
deformation is that due to the scan angle effects, and the different 
scales in the x and y directions. For example, the triangular area 
(arrow A in Figures 4.2 and 4.5) near the top of flight 208 formed 
by the intersecting roads and the stream cutting across them has 
been visibly deformed, showing a widening at the top and a resultant 
lessening of the inclination with which the stream cuts across the 
flight line in this area. Another good example of these deformations 
is found near the bottom of flight 208. Here an interstate highway 
cuts diagonally across the flight line (arrow B, Figure 4.2). 

Careful visual study reveals image deformation displacing the originally 
straight alignment of the roadway into a curved alignment. Also 
apparent, particularly on the display for flight 208 (Figure 4.2), 
is the perturbation due to changes in Y c , the sensor Y position. 
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as evidenced by the departure of the central roadway from the straight 
alignment shown on the map sheet. Upon close inspection of the 
display for flight 208 (Figure 4.2), the effects of scan time and/or 
k rotation down the flight line, is evidenced by the rotation of 
the cross roads and field boundaries. Scan angle deformation is 
also apparent in flight 218, near the river bend area. The river 
bend itself shows a visible spread due to these effects (arrow A, 
Figures 4.4 and 4.6). Also the road running diagonally across the 
flight line above the river upon close inspection shows a reduction 
in the inclination (arrow B, Figures 4.4 and 4.6). 

4.2 Flow Chart for MSS Digital Data Analysis 
Figure 4.7 illustrates the system for analysis of digital 
MBS data arrays which was used in this investigation. The flow 
chart shown displays a method of integrating the principles and 
mathematical techniques discussed in the previous chapter into a 
coherent and organized system for the analysis and restitution of 
digitally recorded multispectral scanner arrays. Referring to 
Figure 4.7* the first step in the procedure, after the desired control 
points have been chosen, is the digitization of map data for the 
control point coordinates, and contour information for subsequent 
assignment of element elevations. For this investigation, the 
digitization was carried out on a LARR-V cartesian coordinate 
digitizer at Purdue University. In the following, the coordinate 
axes of the digitizing equipment will be denoted U-V to differentiate 
them from the object space X-Y system, or the image x-y system. 
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In Figure h.Q is shown the operation of this digitization equipment. 
The material from which U-V coordinate values are to be digitized 
(in this case, the map sheet) is fastened upon a large flat bed. 
Points are located by the operator with a cursor operating along 
perpendicular U-V axes. The U-V coordinates under the cursor are 
instantaneously displayed upon a console shown in Figure if. 8 directly 
behind the operator. The U-V coordinate values of a point may be 
automatically recorded upon punched cards using the keypunch shown 
in the lower right foreground of Figure if. 8. This keypunch is 
interfaced with the digitizer, and the operator activates a relay 
when he is over a point for which he wants the U-V coordinate value 
recorded on cards. 

For this investigation, the following steps were followed in 
the digitization of the map information. 

1. ) The map sheet was fastened upon the flat bed. An attempt 

was made to approximately align the direction of the flight 
line with the U axis of the digitizer. 

2. ) Universal Transverse Mercator (UTM) grid marks on the 

map edges were digitized, first along the map edge aligning 
more nearly to the digitizer U axis, then along the other 
edge. These 1000 m (3281 ft.) intervals could then be 
used to define an absolute scale to transform data arrays 
to ground system coordinates after the analysis and/or 
restitution has been accomplished. 

3. ) All reference point U-V coordinates were digitized in the 


order of their numbering. 
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4.) Contour information was digitized. The method used was 
the contour following mode, in which selected contours 
were followed and U-V coordinates were digitized at random 
sample intervals along each contour. These values were 
recorded directly upon punched cards. Each card was formatted 
to record first a count index, representing the total 
number of contour points recorded previously for the contour 
being followed. Next, the elevation of the contour was 
recorded, followed by 8 values of U-V coordinates of 8 
of the sample points digitized along the contour. Using 
this format, a total of about 1000 cards resulted for 
flight 208, and approximately 3000 cards were required 
for the more rugged terrain of flight 218. 

After digitization of the map data, this information was trans- 
formed from the U-V digitizer system into a cartesian coordinate 
system at the image, or array scale. The transformation used to 
generate these scaled control and contour X-Y points was the well 
known conformal similarity transformation. The equations of the 
transformation are of the form 

X - X(U cos <j> + V sin $) + X 

(4,1) 

Y = X( -U sin 4> + V cos <{>) + Y c 

in which Xq, Y 0 are translations, <J> is the rotation between the 
two systems, and X is a scale factor. These four parameters Were 
determined using least squares techniques, to bring the digitized 
map data to the same approximate scale and orientation as the digital 
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data arrays. All coordinates of the digitized map data were then 
transformed using Equation (4.1), and all subsequent steps were 
performed at image scale. Table 4.1 shows the values of the reference 
point positions for flight 208 after their transformation to image 
scale, and the line-column array positions for the corresponding 
points within the data array. Table 4.2 shows the reference point 
coordinates for flight 218. After these preliminary steps, the 
elements of the arrays, were transformed to the scaled control values 
using an affine transformation, and elevations were assigned to each 
element. Geometric analysis and restitution was then performed 
using the methods described in Sections 3.1 - 3.3. These steps 
in the flow chart are explained in more detail in the following 
sections. 

4.3 Results of Elevation Assignment Algorithms 

In Section 3.4, several possible methods of assigning elevations 
to elements within the data arrays were presented. Two of the 
possible methods were tested in this investigation. For flight 208, 
having flat terrain and hence relatively few digitized contour 
points, both the contour-scan and arithmetic mean nonpararaetric 
methods were attempted. For flight 218, having more rugged terrain 
and many more digitized contour points, only the contour-scan inter- 
section method was done, for reasons which will be discussed sub- 
sequently. 

As pointed out in Section 3.4.1, the affine transformation was 
used within the element elevation assignment program, to approximately 
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Table 4.1 



Reference 

Point Coordinates 

for Flight 208 




Control Points 



Point 

Transformed 



Number 

Map Position 

Array Position 


X 

Y 

Line 

Col 




(x) 

(y) 

3 

209.3 

47.1 

215 

26 

7 

210.8 

164.6 

208 

187 

8 

352.9 

47.0 

356 

29 

11 

442.1 

47.6 

447 

27 

13 

559.6 

47.9 

564 

26 

15 

593.9 

105.9 

595 

104 

17 

677.4 

106.0 

679 

103 

18 

612.3 

194.2 

611 

218 

20 

793.6 

106.7 

798 

105 

22 

903.7 

49.6 

909 

27 

24 

908.8 

164.9 

907 

1 86 

25 

964.6 

106.7 

965 

106 

28 

1024.1 

164.2 

1024 

185 

29 

1136.1 

51.4 

1139 

31 

30 

n4i.i 

164.9 

1138 

187 

31 

1253.6 

50.8 

1254 

30 

34 

1316.1 

166.7 

1310 

189 

38 

1490.6 

167.1 

1479 

194 

4o 

34.7 

104.6 

30 

97 

43 

57.9 

192.6 

45 

215 

44 

91.8 

46.9 

93 

21 

47 

152.2 

104.9 

153 

105 

54 

325.4 

104.8 

324 

107 

55 

325.7 

164.0 

321 

188 

63 

443.3 

104.6 

447 

104 

64 

444.0 

193.9 

443 

218 

66 

499.0 

105.2 

502 

103 

68 

502.3 

193.6 

501 

218 

71 

624.1 

46.7 

630 

26 

72 

727.7 

164.9 

726 

186 

74 

734.5 

47.4 

741 

27 

75 

846.8 

48.6 

854 

27 

77 

850.9 

164.3 

852 

187 

83 

1048.5 

106.4 

1051 

105 

89 

1308.3 

108.3 

1305 

107 

92 

1428.6 

21.1 

1426 

1 

94 

1488.9 

110.1 

1482 

116 

97 

1571.7 

48.0 

1568 

33 

100 

1571.4 

166.8 

1565 

195 
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Point 

Number 


1 

2 

4 

5 

6 

9 

10 
12 
14 
1 6 
19 
21 
23 
26 
27 

32 

33 

35 

36 

37 
39 

41 

42 

45 

46 

48 

49 

50 

51 

52 

53 

56 

57 

58 

59 

60 

61 

62 
65 
67 
69 


Table 4.1, cont. 

Check Points 

Transformed 

Map Position Array Position 


X 

Y 

Line 

Col 



(x) 

(y) 

92.5 

105.2 

90 

102 

179.7 

47.1 

185 

26 

209.7 

75.0 

212 

62 

210.5 

105.4 

212 

105 

92.9 

132.6 

88 

l4l 

353.7 

69.1 

355 

64 

354.5 

105.3 

355 

106 

443.6 

165.2 

443 

186 

560.4 

105.8 

564 

104 

677.1 

47.3 

681 

26 

792.4 

48.4 

801 

28 

795.2 

164.2 

795 

185 

905.8 

106.6 

908 

105 

1022.0 

106.4 

1024 

106 

1023.3 

136.0 

1024 

148 

1311.9 

50.5 

1310 

30 

1314. 0 

108.3 

1309 

107 

1370.7 

50.2 

1370 

27 

1486.9 

50.9 

1484 

33 

1492.6 

50.4 

1490 

33 

1490.2 

138.3 

1482 

157 

35.9 

162.9 

28 

179 

63.5 

75.5 

62 

58 

122.2 

77.0 

122 

64 

123.0 

104.8 

121 

105 

152.6 

134.0 

149 

148 

180.5 

74.9 

184 

62 

181.3 

104.9 

182 

106 

210.4 

134.2 

210 

148 

268.7 

164.3 

263 

189 

325.0 

46.9 

326 

29 

384.9 

164.6 

384 

188 

404.6 

46.7 

410 

27 

397.2 

164.6 

396 

187 

412.9 

47*1 

4l8 

27 

412.9 

75.0 

4i8 

62 

414.0 

164.7 

415 

187 

442.5 

75.0 

447 

60 

453.9 

104.7 

457 

104 

500.7 

134.8 

503 

146 

530.8 

105.3 

535 

104 
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Table 4.1, cont. 

Check Points 

Point Transformed 

Number Map Position Array Position 

X Y Line Col 

(x) (y) 


73 

764.5 

105.4 

764 

104 

76 

848.8 

106.0 

851 

105 

78 

879-2 

134.8 

880 

148 

79 

879-6 

164.4 

879 

186 

80 

936.3 

106.2 

936 

106 

81 

1020.3 

78.5 

1025 

64 

82 

1029.6 

78.6 

1032 

65 

84 

1110.7 

135.7 

1110 

149 

85 

1137-8 

77-9 

1138 

65 

86 

1194.9 

49.8 

1198 

29 

87 

1197. ^ 

107.6 

1198 

107 

88 

1256.1 

112.6 

1253 

117 

90 

1315-3 

136.6 

1310 

150 

91 

1344.4 

136.7 

1339 

150 

93 

1429-4 

49.9 

1425 

31 

95 

1545-5 

78.8 

1541 

71 

96 

1546.4 

110.2 

1538 

120 

98 

1571-6 

78.9 

1566 

74 

99 

1571.5 

110.2 

1566 

120 

Note: 

Check point 70 was 

omitted after 

a mistake 

in digitization 


for this point was 

discovered. 
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Table 1+ .2 

Reference Point Coordinates for Flight 218 
Control Points 


Point Transformed 

Number Map Position Array Position 



X 

Y 

Line 

Col 




(x) 

(y) 

1 

5.1+ 

88.5 

9 

91 

3 

115.7 

166.7 

117 

178 

1+ 

132.6 

86.3 

137 

86 

6 

269.5 

121+.0 

265 

130 

7 

313.7 

138.7 

310 

1 1+6 

9 

; 1+1+1. 3 

51.6 

1+1+5 

36 

11 

1+83.8 

183.8 

1+85 

190 

12 

1+92.9 

127.8 

1+95 

128 

13 

529.3 

93.8 

532 

85 

11+ 

588.1 

37.2 

589 

26 

16 

690.0 

1+2.2 

695 

26 

17 

688.0 

103.7 

692 

100 

19 

871.2 

121+.2 

877 

12 1+ 

20 

928.1 

212.1 

922 

219 

22 

993.3 

130.1 

99 k 

136 

23 

999.5 

61.5 

1001 

52 

2k 

1101.7 

122.9 

1099 

131 

2; 

1116.1 

16.2 

1122 

12 

28 

1171+.9 

126.7 

1171 

131+ 

29 

1321+.6 

33.8 

1325 

24 

30 

1296.5 

206.2 

1289 

213 

32 

1I+07.I+ 

66.6 

11+09 

57 

33 

1362.0 

177.7 

1358 

189 



Check Points 



2 

66.1 

113.0 

63 

117 

5 

257.6 

108.1 

256 

109 

8 

1+36.8 

172.8 

1+37 

182 

10 

1+1+6. 6 

76.5 

1+1+9 

65 

15 

627.5 

99.1+ 

627 

9l+ 

18 

71+7.7 

120.7 

7I+8 

121 

21 

985.8 

133.9 

986 

ll+2 

27 

1121 . 3 

106.8 

1120 

no 

31 

1367.5 

61+.5 

1372 

56 


Note: Control point 26 was omitted after a mistake in digitization 

for this point was discovered. 
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transform the data arrays into the conformal system of the scaled 
control points. Parameters for this transformation were determined 
by least squares methods utilizing as input the scaled control point 
coordinates and the data array positions for the control points, 
as given in Tables 4.1 and 4.2. The flight lines were divided into 
sections, and an independent determination of the affine parameters 
was determined within each flight line section for the affine equa- 
tions of the form 


X » A 0 + A-j_ x + Ag y 
Y - B Q + x + B2 y 


(4.2) 


Table 4.3 portrays the results of this least squares fit for flight 
208. In this instance, the flight line was divided into five sections. 
Table 4.4 depicts the results of the analysis for flight 218, in 
which the flight line was segmented into four sections . The resulting 
standard deviations and variances are given in terms of array 
elements. The average ground size of such an element for the arrays 
tested is about 40 ft. diameter. This unit is often called a remote 
sensing unit, and is abbreviated rsu. Therefore, the statistical 
quantities given in such units reflect the accuracies of analysis 
methods in terms of the resolving power of the system, rather than 
absolute distances. 

The pooled variance in these tables was computed as the summation 
for all sections divided by the summation of the degrees of freedom 


for all sections. 



Table 4.3 


Statistics from Affine Tran sformat ion for Flight 208 


Section Number 

1 

2 

3 

4 

5 

Pooled* 

Degrees of Freedom 

l4 

Ik 

12 

12 

16 

68 

A posteriori Estimate of 
Variance ( rsu 2 ) 

3.71 

2.1k 

2.84 

0.38 

H 

00 

$ 

M 

2.20 

A posteriori Estimate of 
Standard Deviation (rsu) 

1.93 

1.46 

1.69 

0.6l 

1.34 

1.48 

Statistics 

Section Number 

from Affine 
1 

Table 1+.1+ 
Transformation for 
2 3 

Flight 218 
4 

Pooled* 


Degrees of Freedom 

8 

12 

12 

16 

48 


A posteriori Estimate of 
Variance (rsu 2 ) 

5.91 

7.18 

10.35 

12.26 

9.45 


A posteriori Estimate of 
Standard Deviation (rsu) 

2.1+3 

2.68 

3.22 

3.50 

3.07 



k 

E (d.o.f.^ (t^ 2 ) 

* pooled o 2 » i=l = where k is number of sections 

k 

I (d.o.fJi 
i=l 



The parameters thus generated were used within the element 
elevation assignment program to approximately fit the data arrays 
to the scaled control points. To test the suitability of this 
transformation, the analysis described in Section 3.4.1 must be per- 
formed. Ground slopes within flight 208 were very small, rarely 
exceeding 0.1, and the resultant allowable error for the affine 
transformation from Equation (3.53) is on the order of 12 resolution el- 
ements, even near the strip edges where the total scan angle approaches 
38.5°. Clearly, the a posteriori reference standard deviations 
obtained from the adjustment yield values well within this limit. 

For flight line 218, ground slopes were much steeper, approaching 
0.5 in some places within the flight line. Assuming the worst case, 
in which the maximum slope of 0.5 occurs at the maximum scan angle 
of 38.5° » the allowable error is 2.5 resolution elements. The 
resultant a posteriori pooled reference standard deviation for this 
flight line was 3.07 resolution elements. This means that after element 
elevations were assigned, and subsequent analysis was done, some 
of the residual error within flight line 218 would possibly be 
due to topographic relief effects because of errors in assignment 
of element elevations. These values however, represent the worst 
case . Over most of the flight line , allowable planimetric error 
would be larger than the 2.5 resolution elements given, and resulting 
relief displacement errors would be reduced. For scan angles less 
than 26.5° for example, the allowable planimetric error would be 
greater than 4 resolution elements * 
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In Figures 4.9 and 4.10 are shown typical profiles from each 
of the two flight lines from which an indication of the reliability 
of the height assignment algorithms may be gained. In Figure 4.9 
is depicted a typical profile across the flight line for data array 
line 7^3 of flight 208. In this figure, the solid line represents 
a reference ground profile taken from the map sheet. The circle- 
dash combination represents the profile obtained by the contour- 
scan intersection method of element elevation assignment, and the 
dashed line represents the profile resulting from element elevation 
assignment using the arithmetic mean method. It should be noted 
that the elevations from the map are also subject to some error. 

The National Map Accuracy Standards for vertical map information 
state that 90 % of the elevations interpolated from map contours 
shall be correct within 1/2 contour interval. For the relatively 
flat terrain in this flight line, the contour interval was 5 ft. 

( 1.52 m) resulting in a 90 % error of 2.5 ft. (O .76 m) and an approximate 
standard deviation of only 1.5 ft. (0.46 m) due to the map. It may- 
be seen in Figure 4.9 that the profile from the map in almost all 
cases differs from that obtained by either of the elevation assignment 
methods by less than ten feet (3.05 m) . In Figure 4.10, a similar 
portrayal of a typical scan line for flight 218 is shown. For this 
flight line, the map contour interval was 10 ft. ( 3.05 m) and elevation 
errors for points taken from the map were therefore 5 ft. ( 1.52 m) 
at the 90 $ level, for an approximate standard deviation of 3 ft. 

(0.92 m). Mismatches of up to 40 ft. (12.2 m) between the map profile 


Elevation (ft.) 
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and that obtained by element elevation assignment are seen to exist 
for short sections along the scan line. As will be shown subsequently, 
the allowable elevation assignment error for these flight lines 
is on the order of 60 feet. Thus, the maximum elevation assignment 
errors for these sample profiles (10 feet and 40 feet for flights 
208 and 218, respectively) would result in planimetric errors of less 
than one resolution element . 

A statistical test was performed to evaluate the sufficiency 
of the element elevations assigned. For this test, the elevations 
assigned to the reference point array locations by the element 
elevation assignment programs were compared with elevations for the 
corresponding points taken from the map sheet. The differences 
between elevations assigned from the programs and those taken from 
the map sheet were used as the data upon which to perform the statis- 
tical analyses. Table 4.5 gives the data generated for the control 
points of flight 208. In each case, the column noted as ’'Difference’' 
shows the value obtained by subtracting the value taken from the 
map sheet from that assigned by the elevation assignment algorithm 
noted. Table 4.6 gives the same information for flight 218. 

In order to evaluate the results of a statistical analysis, 
it is important to state the objectives and the desired results 
of such an analysis. First, it would be desired that the mean of 
the elevation differences be statistically near zero, in order that 
it may be stated that, on the average, the height assignment al- 
gorithms are assigning correct elevations. In addition it is desired 
that the variance, and the resulting standard deviation, be such 


Table 4.5 

Elevation Information for Flight 208 



A. 

B. 

C. 

D. 

E. 

Control 

Elevation 

Elevation Assigned 


Elevation 


Point 

from Map, 

by Contour-Scan, 

Difference , 

Assigned by 

Difference, 

Number 

in ft. (m) 

in ft. (m) 

B-A in ft. (m) 

Arithmetic 

D-A in ft. (m) 





Mean, in ft. 

(m) 

3 

695 (211.8) 

695 (211.8) 

0 (0) 

696 (212.1) 

1 (0.3) 

7 

705 (214.9) 

705 (214.9) 

0 (0) 

705 (214.9) 

0 (0) 

8 

709 (216.1) 

710 (216.4) 

1 (0.3) 

710 (216.4) 

1 (0.3) 

11 

708 (215.8)* 

710 (216.4) 

-2 (0.6) 

710 (216.4) 

2 (0.6) 

13 

706 (215.2) 

709 (216.1) 

3 (0.9) 

710 (216.4) 

4 (1.2) 

15 

685 (208.8) 

695 (211.8) 

10 (3.0) 

696 (212.1) 

11 (3.3) 

17 

704 (214.6) 

705 (214.9) 

1 (0.3) 

704 (214.6) 

0 (0) g 

18 

719 (219.2) 

720 (219.5) 

1 (0.3) 

715 (217.9) 

-4 (-1.2) 

20 

715 (217.9)* 

720 (219.5) 

5 (1.5) 

719 (219.2) 

^ (1.2) 

22 

723 (220.4) 

720 (219.5) 

-3 (-0.9) 

717 (218.5) 

-6 (-1.8) 

24 

715 (217.9)* 

715 (217.9) 

0 (0) 

715 (217.9) 

0 (0) 

25 

718 (218.8)* 

715 (217.9) 

-3 (-0.9) 

715 (217.9) 

-3 (-0.9) 

28 

722 (220.1) 

720 (219.5) 

-2 (-0.6) 

718 (218.8) 

-4 (-1.2) 

29 

715 (217.9) 

715 (217.9) 

0 (0) 

715 (217.9) 

0 (0) 

30 

710 (216.4)* 

720 (219.5) 

10 (3.0) 

710 (216.4) 

0 (0) 

31 

721 (219.8) 

724 (220.7) 

3 (-0.9) 

722 (220.1) 

1 (0.3) 

34 

731 (222.8) 

720 (219.5) 

-11 (-3.3) 

721 (219.8) 

-10 (3.0) 

38 

700 (213.4)* 

711 (216.7) 

11 (3.3) 

713 (217.3) 

3 (-0.9) 

40 

751 (228.9) 

755 (230.1) 

4 (1.2) 

751 (228.9) 

0 (0) 

43 

751 (228.9) 

746 (227.4) 

-5 (1.5) 

753 (229.5) 

2 (0.6) 

44 

741 (225.9) 

745 (227.1) 

4 (1.2) 

731 (222.8) 

-10 (3.0) 

47 

743 (226.5) 

740 (225.6) 

-3 (-0.9) 

744 (226.8) 

1 (0.3) 

54 

745 (227.1) 

745 (227.1) 

0 (0) 

744 (226.8) 

-1 (-0.3) 

55 

741 (225.9) 

740 (225.6) 

-1 (-0.3) 

743 (226.5) 

2 (0.6) 

63 

733 (223.4) 

731 (222.8) 

-2 (-0.6) 

732 (223.1) 

-1 (-0.3) 

64 

729 (222.2)* 

730 (222.5) 

-2 (-0.6) 

727 (221.6) 

-2 (-0.6) 



Table 4.5, cont. 



A. 

B. 

C. 

D. 

E. 

Control 

Elevation 

Elevation Assigned 


Elevation 


Point 

from Map, 

by Contour-Scan, 

Difference, 

Assigned by 

Difference , 

Number 

in ft. (m) 

in ft. (m) 

B-A in ft. (m) 

Arithmetic 
Mean, in ft. (m) 

D-A in ft. (m) 

66 

730 (222.5) 

730 (222.5) 

0 (0) 

730 (222.5) 

0 (0) 

68 

725 (221.0) 

730 (222.5) 

5 (1.5) 

725 (221.0) 

0 (0) 

71 

731 (222.8) 

730 (222.5) 

-1 (-0.3) 

730 (222.5) 

-1 (-0.3) 

72 

723 (220.4) 

730 (222.5) 

7 (2.1) 

726 (221.3) 

3 (-0.9) 

74 

731 (222.8) 

725 (221.0) 

-5 (-1.5) 

730 (222.5) 

-1 (-0.3) 

75 

729 (222.2) 

725 (221.0) 

-4 (-1.2) 

728 (221.9) 

-1 (-0.3) 

77 

727 (221.6) 

725 (221.0) 

-2 (-0.6) 

725 (221.0) 

-2 (-0.6) 

83 

730 (222.5)* 

725 (221.0) 

-5 (-1.5) 

730 (222.5) 

0 (0) 

89 

732 (223.1) 

730 (222.5) 

-2 (-0.6) 

730 (222.5) 

2 (0.4) 

92 

740 ( 225 . 6 )* 

740 (225.6) 

0 (0) 

740 (225.6) 

0 (0) 

94 

741 (225.9) 

740 (225.6) 

-1 (-0.3) 

740 ( 225 . 6 ) 

1 (0.3) 

97 

729 (222.2) 

735 (224.0) 

6 (1.8) 

731 (222.8) 

2 (0.6) 

100 

731 (222.8)* 

736 (224.3) 

5 (1.5) 

731 (222.8) 

0 (0) 


* These values were noted on the map as spot elevations. 

Other values were obtained by interpolating between contours 
on the map. 


Resulting Statistics 

Data Set C Data Set E 

Degrees of freedom 38 38 

Mean, in ft. (m) +0.56 (+.17) -0.15 (-0.4) 


Variance 


20.73 


12.77 
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Table 4.6 


Elevation Information for Flight 218 



A. 

B. 

c. 

Control 

Elevation 

Elevation Assigned 


Point 

from Map, 

by Contour -Sc an , 

Difference 

Number 

in ft. (m) 

in ft. (m) 

B-A in ft. 

1 

595 (181.4) 

599 (182.6) 

4 (1.2) 

3 

559 (170.4) 

625 (190.5) 

66 (20.1) 

4 

587 (178.9) 

609 (185.6) 

22 (6.7) 

6 

615 (187.4) 

599 (182.6) 

-16 (-4.9) 

7 

580 (176.8) 

599 (182.6) 

19 (5.8) 

9 

705 (214.9)* 

699 (213.1) 

-6 (-1.8) 

11 

561 (171.0)* 

574 (175.0) 

13 (4.0) 

12 

587 (178.9) 

599 (182.6) 

12 (3.7) 

13 

673 (205.1)* 

663 (202.1) 

-10 (-3.0) 

14 

593 (180.7)* 

6l4 (187.1) 

21 (6.4) 

16 

715 (217.9) 

690 (210.3) 

-25 (-7.6) 

17 

725 (221.0) 

687 (209.4) 

-38 (-11.6) 

19 

759 (231.3) 

748 (228.0) 

-11 (-3.4) 

20 

660 (201.2) 

649 (197.8) 

-11 (-3.4) 

22 

650 (198.1) 

649 (197.8) 

-1 (-0.3) 

23 

581 (177.1) 

599 (182.6) 

18 (5.5) 

24 

580 (176.8) 

574 (174.9) 

-6 (-1.8) 

25 

543 (165.5) 

550 (167.6) 

7 (2.1) 

28 

633 (192.9) 

625 (190.5) 

-8 (-2.4) 

29 

725 (221.0) 

717 (218.5) 

-8 (-2.4) 

30 

632 (199-6) 

550 (167.6) 

-82 (-25.0) 

32 

552 (168.2) 

576 (175.6) 

24 (7.3) 

33 

613 (186.8)* 

549 (167.3) 

-64 (-19.5) 


(m) 


* These values -were noted on the map as spot elevations. Other 
values were obtained by interpolating between contours on the map. 


Resulting Statistics 


Degrees of freedom 
Mean, in ft. (m) 
Variance 


Data Set C 
22 

-3.48 (- 1 . 06 ) 
924.1 
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that the great majority of the data array positions would be within 
the allowable height assignment error as derived in Section 3.4.1. 
Referring to Figure 3*5 and Equation (3.51) > it may be seen that 
the minimum value of this allowable error would occur near the flight 
line edge, at a scan angle of 38.6°. At this scan angle, with a 
flying height of 5000 ft. (1.52 km), the ground size of a resolution 
element in the Y direction is approximately 50 ft. (l4.9 m) and the 
resulting allowable height assignment error is approximately 6l ft. 
(18.8 m) . At a scan angle of 20°, the resulting ground size of a 
resolution element is 34 ft. (10.4 m) and the resulting allowable 
height assignment error is 93 ft. (28.5 m) . Thus, it may be seen 
that 60 ft. would represent a conservative value, the square of which 
yields a variance against which to test the results from the height 
assignment algorithms. 

A complete statistical analysis of flight 208 was not performed. 

A cursory inspection of the differences between the assigned element 
elevations and those from the map sheet (data sets C and E, Table 4.5) 
reveals that no value in the data even approaches the minimum accep- 
table value of 60 feet. Hence for this flight line the elevations 
may be assumed to be assigned with suitable accuracy, using either 
the contour-scan or arithmetic mean method. Due to its computational 
efficiency, the contour-scan method was utilized in subsequent analysis 
and restitution computations when element elevations were required. 

For flight 218, the data from which statistical quantities 
were computed for testing were the differences between the elevations 
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assigned and those taken from the map (data set C, Table 4.6). 

As a first step, the data were tested for normality, in order that 
subsequent tests upon the mean and the variance would be valid. 

The Shapiro-Wilk test for small samples was used [108], In this 
test a statistic, denoted w, is computed and compared with tabular 
values published in that reference in order to test the hypothesis 
that the data are normally distributed against the alternative 
that the data are not normally distributed. When the computed w 
from the data is below the appropriate tabulated value, then non- 
normality of the data is indicated. The value of w calculated from 
the data of flight 218 was 0.925. The value of w is computed as 

_b£_ 
w = s.s. 

in which s.s. represents the sum of squares of the data values. 

The quantity b is calculated by first ordering the data in ascending 
order (d ls d2» . . . dn where is the smallest value and the 
largest data value). If n is an even integer, then a parameter k 
is defined such that n * 2k and b is computed with 

k 

b = 2 ®n-i+i (^n-i+l ** ^ 

i-1 

where the a values are coefficients tabulated in the reference. 

If n is an odd integer then the central data value is omitted from 
the set and k is defined such that n = 2k + 1 and the expression 
for b becomes 

b = a n (d n - d x ) + . . . + aj^dj^g - d k ) 
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where the value of d^+i, the sample median, does not enter the com- 
putation of b. The tabular values for the appropriate data sample 
size (23) is 0.91^ for an a level of 0.05 and 0.928 for an a level 
of 0.10. The hypothesis that the data are normal may then be accepted 
at an a level of 0.0?. 

Next, the mean value of these data was computed and tested. 

The null hypothesis formulated was that the mean was equal to zero, 
versus the alternate hypothesis of a non-zero mean. The hypothesis 
was tested at an a level of 0.05 • At statistic was calculated 
from the data and found to be -0.5^9. The tabulated value for 
■t.975 22 was to be 2.074. Thus, for a two-sided test, the 

null hypothesis would be rejected only if the t value calculated 
was greater than 2.07^ or less than -2.07^4, and the hypothesis of 
zero mean was accepted. 

The final test was performed upon the variance. This represents 
an important step, since it is here that the adequacy of the height 
assignment is determined. For this step, the sample variance of 
92^.1 calculated from the data was compared with the allowable height 
assignment variance assumed as 3600 (from the allowable height as- 
signment error of 60 feet) . It was assumed that the degrees of freedom 
associated with the allowable height assignment error is infinite. 

The null hypothesis was formulated that the variance calculated from 
the data was less than or equal to 3600, versus the alternative 
that this variance is larger than 3600. The hypothesis may then 
be tested using a x 2 statistic, or an F statistic having infinite 
degrees of freedom in the denominator. A x 2 statistic was calculated 


4 ->fV; 
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from the data and found to be 5.65. The test was performed at an a 
level of 0 . 05 , for which the tabulated value of x 2 95 22 was ^ oun ^ 
to be 33.9. The null hypothesis would be rejected only if the x 2 
value calculated from the data was greater than the value tabulated, 
and the hypothesis was accepted. 

The net result of this test is that the elevation data assigned 
to the element array positions may reasonably be assumed to lie 
within 60 ft. of the true value. The resulting planimetric error 
due to elevation assignment errors will then be less than one res- 
olution element. With this knowledge, attention is directed in the 
next section to the problem of testing the various methods of analysis 
and restitution previously presented. 

4.4 Comparison of Analysis Methods 

The purpose of this section is to systematically test the analysis 
and restitution methods presented in Chapter 3, when operating upon 
the actual data arrays chosen. For each method, the assumptions 
made for each analysis will be presented, numerical data resulting 
from the analysis will be given, and meaningful differences in the 
results under different assumptions will be pointed out. 

4.4.1 Collinearity Equations 

In testing this method a two step procedure was carried out 
in which control points from the map and the corresponding array- 
positions were used to estimate parameters of the exterior orientation, 
followed by computation of coordinates of check points withheld 
from the adjustment. This results in two useful statistics. 
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1. ) The a, posteriori reference variance which is computed from 

the adjustment itself. 

2. ) A "positional check variance", based upon the variances 

in the X and Y directions computed from the check points. 

For this purpose the square of the circular standard error as defined 
in [113] was used. A conservative estimate for this value c, is 
given in this reference as c = 0.5(a x + a ). These quantities 

*J 

were assigned n - 1 degrees of freedom, where n is the number of 
check points, and were statistically tested for significant differences, 
in order to ascertain which of the assumed functional forms produced 
better results. 

For each flight line, several cases with the following functional 
assumptions were performed, and the results from each were compared. 

Case Cl: For this case, the functional forms for the exterior 

orientation parameters were of the form 

x c = 

Y c = 

Zc 52 

K ~ 

$ = 

0 ) = 

In which the recoverable orientation elements are assumed as linear 
functions of x within a flight line section. 


A x + B x x 
Ay + By X 
A z + B z x 


A K + B k x 


(4.3) 
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Case C2: For this case, functional forms are assumed to be 


X c - A x + B x x + C x x 2 

Y c = Ay + By x + Cy X 2 

Z c = + B z x + C z x 2 

k - A^, + B K x + C K x 2 

<{> = 0 
u) = 0 


(h.k) 


These polynomials were used to investigate the usefulness of non- 
linear variations of orientation elements within a flight line section 
Case C3: For this case , at. attempt was made to assume polynomial 

forms which would correct for non-linear effects 
within a section, and retain computational efficiency 
for the later element by element transformation of 
the arrays. The functional forms assigned here were: 


X c = A x + B x x + Cy x 2 

Y c = Ay + By X + Cy X 2 


Z c = 
K 

4> 

to 


A z + B z x 


(4.5) 


= A k 
■ 0 

= 0 


Using these polynomials , the intent was to model the non-linear 
effects within each section using non-linear functions which affect 
most directly the element positions, i.e. X c , Y c . The k tern is 
modeled as a constant within each section in order to avoid the 
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re-evaluation of trigonometric functions for each scan line. Note 
here that only a linear form was assumed to model the flying altitude. 
Experience gained during the investigation showed that the addition 
of an x 2 term in the polynomial used to model this element produced 
negligible improvement over the case using the functional form shown. 

Case C4: In this case the following polynomial forms were 

assumed: 


X c = Ay + By x + Cy x 2 

Y c = Ay + By x + Cy x 2 

Z c = Ag + B Z * 

K = Ay, + B K X 
4 ) = 0 

oj = 0 


(4.6) 


The purpose here was to test if a significant difference could 
he detected by linear modeling of the k term over the computationally 
more efficient model in which tc is assumed constant for the flight 
line ( Case C3) • 

For each of the cases given above, three different procedures 
were used. The adjustment and checking was done first treating 
the entire flight line as a unit. That is, it was assumed that a 
single set of parameters was valid for the entire length of the 
flight line. Then the flight line was treated in two sections, 
in which a separate set of parameters was solved for in each section. 
For this mode, the flight line was divided into sections each about 
half the length of the flight line. As a third case, the flight 
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line was segmented into thirds and separate parameter sets computed 
for each of the three sections. For the instances in which the flight 
line was divided into sections, the constraints of Equation set (3.39) 
were enforced at the section boundaries. 

A summary of the results for flight 208 is shown in Table 4.7* 

For each case, the total number of parameters is given, and the 
number of degrees of freedom for the adjustment. The reference 
variance resulting from the adjustment is next shown. This reference 
variance is a unitless quantity which, when multiplied by the assumed 
a, priori variances of the observations will result in a posteriori 
estimates of these variances, as indicated by the least squares 
adjustment. For flight line 208, the a priori estimates of variance 
assumed were 1.0 for the scaled control points from the map, and 2.25 
for the corresponding array positions taken from the line printer 
display. Table 4.8 shows the corresponding information for flight 218. 
For this flight line, the a priori estimates of variance on the 
control points were 1.0 for the scaled X-Y map coordinates and 6.25 
for the corresponding x-y array positions. 

The units for the standard deviations are in terms of elements 
within the data arrays, or remote sensing units (rsu). The estimates 
of standard deviation taken for the digitized map control points 
were based upon the National Map Accuracy Standards. The standards 
state that 90 % of all "well defined" points plotted fall within 
1/50 in. (0.5 mm) of their true position. Assuming equal errors 
in X and Y, the component 90% error estimates would be E^ = E^ * E/»^" 
where E^, Ey are the 90% errors in X and Y, and E is the positional 


226 


Table 4.7 

Statistics from Collinearity Analysis of Flight 208 
(60 Check Points) 


Number of Sections 


12 3 


Case Cl: Number of Parameters 


8 

16 

24 

Degrees of Freedom 


70 

66 

62 

Reference Variance 


2.58 

l.l4 

0.96 

Check Variance x in rsu 2 (d.o.f. = 

59) 

9.10 

3.05 

2.40 

Check Variance y in rsu 2 (d.o.f. = 

59) 

4.55 

4.04 

3. 4i 

Positional Check Variance in rsu 2 


6 .63 

3.53 

2.88 


Case C2: Number of Parameters 


12 

24 

36 

Degrees of Freedom 


66 

58 

50 

Reference Variance 


1.10 

0.78 

0.59 

Check Variance x in rsu 2 (d.o.f. - 

59) 

2.48 

2.27 

2.03 

Check Variance y in rsu 2 (d.o.f, = 

59) 

4.13 

2.81 

1.86 

Positional Check Variance in rsu 2 


3.25 

2.53 

1.94 


Case C3: Number of Parameters 


9 

18 

27 

Degrees of Freedom 


69 

64 

59 

Reference Variance 


1. 14 

0.86 

O.65 

Check Variance x in rsu 2 (d.o.f. = 

59) 

2.63 

2.53 

2.21 

Check Variance y in rsu 2 (d.o.f. = 

59) 

4.05 

2.81 

1.90 

Positional Check Variance in rsu 2 


3.29 

2.67 

2.06 


Case C4: Number of Parameters 


10 

20 

30 

Degrees of Freedom 


68 

62 

56 

Reference Variance 


1.10 

0.81 

0.57 

Check Variance x in rsu 2 (d.o.f. “ 

59) 

2.56 

2.45 

2.08 

Check Variance y in rsu 2 (d.o.f. = 

59) 

4.06 

2.82 

1.91 

Positional Check Variance in rsu 2 


3.27 

2.64 

1-99 
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Table 4.8 

Statistics from Collinearity Analysis of Flight 218 
(9 Check Points Only) 


Number of Sections 


12 3 


Case Cl: Number of Parameters 


8 

16 

2k 

Degrees of Freedom 


38 

34 

30 

Reference Variance 


1.39 

0.98 

0.87 

Check Variance x in rsu 2 (d.o.f. » 

8) 

7.30 

6.59 

7.30 

Check Variance y in rsu 2 (d.o.f. * 

8) 

4.80 

3.31 

1.89 

Positional Check Variance in rsu 2 


5.98 

4.81 

4.15 


Case C2: Number of Parameters 



12 

24 

36 

Degrees of Freedom 



34 

26 

18 

Reference Variance 



1.08 

0.70 

0.48 

Check Variance x in rsu 2 (d.o. 

f. = 

8) 

6.69 

15.51 

17.86 

Check Variance y in rsu 2 (d.o. 

f. = 

8) 

3.64 

1.63 

0.94 

Positional Check Variance in 

rsu 2 


5.04 

6.80 

6.75 


Case C3: Number of Parameters 


9 

18 

27 

Degrees of Freedom 


37 

32 

27 

Reference Variance 


1.04 

0.79 

0.71 

Check Variance x in rsu 2 (d.o.f. = 

8) 

7.55 

9.92 

8.69 

Check Variance y in rsu 2 (d.o.f. = 

8) 

3.78 

1.91 

1.35 

Positional Check Variance in rsu 2 


5.50 

5.13 

4.22 


Case C4: Number of Parameters 


10 

20 

30 

Degrees of Freedom 


36 

30 

24 

Reference Variance 


1.06 

0.78 

0.65 

Check Variance x in rsu 2 (d.o.f. = 

8) 

6.74 

9.12 

9.44 

Check Variance y in rsu 2 (d.o.f. * 

8) 

3.78 

1.91 

1.34 

Positional Check Variance in rsu 2 


5.15 

4.84 

4.47 
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error given "by the accuracy standards . The resulting standard de- 
viation estimates in X and Y then become approximately .008 in. (0.2 mm) 
each. The resulting ground scale distance when these numbers are 
divided by the 1:24000 scale of the U.S.G.S. J 1/2 minute quadrangle 
maps used is about 10 ft. (3.1 m). This is on the order of one 
half the ground size of a data array element at nadir. However, 
many of the map points digitized were not ’’well defined”, as described 
in the accuracy standards. A value of 1 rsu was therefore assigned 
for the a. priori estimates of standard deviation for map coordinates . 

The a. priori standard deviation for corresponding data array positions 
were based upon the difficulty found in assigning these positions 
within the data arrays. Thus, for flight 208, standard deviation 
estimate of 1.5 rsu (variance 2.25) was chosen, and for the more 
difficult flight 218, a standard deviation of 2.5 rsu (variance 6.25) 
was chosen. 

In order to compare the accuracies obtained for the cases tested, 
tables of F statistics were computed based upon the a posteriori 
reference variances obtained from the adjustments, and the resulting 
positional check variances. These F statistics were computed as 
ratios of resulting variances between the various cases. Table 4.9 
shows the resulting F statistics formed by the quotient of the 
a. posteriori reference variances for the various cases tested for 
flight 208. As an example the value 2.26 in the second row, first 
column of Table 4.9 represents the ratio of the a posteriori reference 
variance of 2.58 obtained in Case Cl treating the entire flight line 
as a unit, divided by the _a posteriori reference variance of 1.14 



Table U.9 


Tabulation of F Statistics for Flight 208 Formed from Reference Variances from Collinearity Analysis 

Case Cl Case C2 Case C3 Case C4 

Number of Sections 



12 3 

12 3 

12 3 

12 3 

Case Cl 1 

2 
3 

1.00 

2 . 26 * 1.00 
2 . 69 * 1.19 1.00 




o 

e 

<D 

O 

ro 

Sections 

W M H 

2 . 34 * 1.04 l.l4 
3.31* 1.46 1.23 
4,37* 1.93* 1.63* 



1.00 

1.4i 1.00 
1 . 86 * 1.32 1.00 



Case C3 b 1 
^ 2 
jg 3 

- - 1 

' 

2.26* 1.00 1.19 
3.00* 1.32 1.12 

3.97* 1.75* 1.48 

1.04 1.46 1 . 93 * 
1.28 1.10 1.46 
1 . 69 * 1.20 1.10 

1.00 

1.32 1.00 

1.75* 1.32 1.00 


a 

Case C4 1 

2 
3 

2.34* 1.04 1.14 
3 . 18 * l.4l 1.18 
4.52* 2 . 00 * 1 . 68 * 

1.00 1.4l 1 . 86 * 
1.36 1.04 1.37 
1.93* 1.37 1.03 

1.04 1.28 1.69* 
1.4l 1.06 1.25 
2.00* 1.51 1.14 

1.00 

1.36 1.00 
1.93* 1.42 1.00 


ro 

ro 

vo 


* Denotes values which are statistically significant at an a level of .05. 








obtained for the same functional assumptions (Case Cl) but dividing 
the flight line into two approximately equal sections. The asterisk 
opposite this value indicates a significant difference between these 
variances. Therefore, significant improvement is noted for this 
data array if linear assumptions are made for the functions of exterior 
orientation and the flight line is divided into two sections over 
the case in which linear functions are assumed and the flight line 
is treated as a single unit. Table 4.10 tabulates the F statistics 
resulting between cases tested if the positional check variances 
are used to form the values. Tables 4.11 and 4.12 show, the resulting 
F statistics obtained in this manner for flight 218. 

Based upon the information in these tables, several facts may 
be pointed out. 

1. ) In general, in order to discern statistically significant 

improvement over the linear assumptions of case Cl, a 
combination of higher order assumptions was necessary, 
coupled with a division of the flight line into three sections. 

2. ) In general, no significant difference was noted between 

Case C2 and Case C3. This suggests that the computationally 
efficient model of Case C3 is adequate for non-linear effects . 
For the arrays tested then, modeling of Z c and k by second 
order functions resulted in no statistical improvement 
over modeling Z c as a linear function and k as a constant 
value . 

3. ) No (statistically significant difference was noted between 

Case C3 and Case c4, suggesting that for the arrays tested, 


Table 4.10 


Tabulation 
Formed from Positional 


of F Statistics 
Check Variances 


for Flight 208 

from Collinearity Analysis 


Case Cl 


Case C2 


Case C3 


Case C4 


Number of Sections 



12 3 

111 

111 

111 

Case Cl 1 

2 
3 

1.00 

1 . 88 * 1.00 
2 . 36 * 1.22 1.00 




tn 

Case C2 ° 1 

-P 2 
0 41 

<D O 
CQ J 

2.04* 1.08 1.13 
2.62 1.39 1.14 
3.41* 1 . 81 * 1.48 

1.00 

1.28 1.00 
1 . 67 * 1.30 1.00 



Vi 

O 

Case C3 1 

u 2 

<t> ^ 

I 3 

2 . 06 * 1.07 1.15 
2.48* 1.32 1.08 
3.23* 1.72* 1.40 

1.02 1.30 1.70* 

1.22 1.05 1.37 
1.59* 1.23 1.06 

1.00 

1.23 1.00 
1 . 60 * 1.30 1.00 


Case C4 1 

2 
3 

2.03* 1.08 1.13 
2 . 52 * 1.34 1.10 
3.32* 1.77* 1-45 

1.00 1.29 1 . 68 * 
1.24 1.04 1.35 
1.63* 1.27 1.03 

1.01 1.22 1 . 59 * 
1.25 1.01 1.28 
1.65* 1.34 1.03 

1.00 

1.25 1.00 
1.64* 1.33 1.00 


* Denotes values vhich are statistically significant at an o level of .05. 
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Table 4.11 


Tabulation of F Statistics for Flight 218 Formed from Reference Variances from Collinearity Analysis 

Case Cl Case C2 Case C3 Case C4 


Number of Sections 



12 3 

12 3 

12 3 

12 3 

Case Cl 1 

2 
3 

1.00 

1.1+2 1.00 

1.60 1.12 1.00 


. 1 ' ■ ' 1 - J - 

1 

m 

Case C2 § 1 

•H _ 

% l 
& 3 

1.29 1.10 1.24 
1.98* 1.1+0 1.21+ 
2.89* 2.01+ i.8l 

1.00 

1.54 1.00 
2.25* 1.46 1.00 



<« 

Case C3 ° 1 

u 2 
<u _ 

1 3 

1.31+ 1.06 1.19 
1.76* 1.21+ i.io 
1.96* 1.38 1.22 

1.04 1,49 2.17 
1.37 1.13 1.64 
1.52 1.01 1.48 

1.00 

1.32 1.00 
1.46 1.11 1.00 


Case C4 1 

2 
3 

1.31 1.08 1.22 
1.78 1.26 1.11 
2.14* 1.51 1.34 



1.02 1.51 2.21* 

1.38 1.11 1.62 

1.66 1.08 1.35 

1.02 1.34 1.49 
1.33 1.01 1.10 
1.60 1.22 1.09 

1.00 

1.36 1.00 
1.63 1.20 1.00 


* Denotes values which are statistically significant at an a level of .05. 
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Table 4.12 


Tabulation of F Statistics 
Formed from Positional Check Variances 

Case Cl Case C2 

Number of 



12 3 

12 3 

Case Cl 1 

2 
3 

1.00 

1.24 1.00 
1.44 1.16 1.00 


0) 

Case C2 q l 

s 2 

& 3 

1.19 1.05 1.22 
1.14 1.4l 1.64 
1.13 1.4o 1.62 

1.00 

1.35 1.00 
1.34 1.01 1.00 

<H 

0 

Case C3 1 

u 2 

1 3 

1.09 1.14 1.32 
1.17 1.07 1-24 
1.42 1.14 1.02 

1.09 1.24 1.23 
1.02 1.32 1.31 

1.20 1.6l 1.60 

Case C4 1 

2 
3 

1.16 1.07 1-24 
1.24 1.01 1.17 
1.34 1.08 1.08 

1.02 1.32 1.31 
1.04 1.40 1.39 
1.13 1.52 1.51 


for Flight 218 

from Collinearity Analysis 


v 

Case C3 Case C4 


Sections 



LO 

u> 
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assuming k to be a constant throughout the flight line is 

practical. 

4. ) For flight 2Q8, having small terrain elevation differences 

and a large number of control points. Tables 4.9 and 4.10 
reveal a statistically significant improvement in the three 
section breakdown as compared with treating the entire 
flight line as a single section regardless of the analytical 
model. This suggests that, for flights having sufficient 
control and low terrain relief differences, adequate seg- 
menting of the flight line should take precedence over the 
use of more complicated analytical models. 

5. ) As shown in Table 4.12, no significant differences between 

the various cases were obtained based upon the check variances 
for flight 218. This was expected since the degrees of 
freedom for these variances was quite small (8). 

6. ) Examination of Table 4.8 points out a danger associated 

with modeling the exterior orientation elements using a 
large number of parameters. For t: is flight line (218) 
having relatively few control points, it may be seen that 
when the higher order assumptions are made to model exterior 
orientation elements, combined with, segmenting the flight 
line, the positional check variance actually increases over 
those cases having fewer parameters. This indicates that 
unless many control points are available, the more highly 
non-linear interpolations may be misleading and dangerous. 
Although these combinations result in lower reference 
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variances from the adjustment, the actual measure of the 
validity of the interpolation, the positional check variance, 
increases, due mainly to the relatively small degrees of 
freedom. 

In summary, the analysis of the data from the two flight lines 
tested suggests that linear modeling of exterior orientation elements 
(Case Cl) coupled with segmentation of the flight line represents 
a modeling method which avoids the problems discussed in 6. above. 
However, the functional forms of Case C3 are recommended. This 
combination of polynomials represents a computationally efficient 
method of accounting for non-linear behaviour in the data, and suitably 
avoids the problems of overparameterization discussed in 5* above. 

This case also results in less computation than for Case Cl, since 
the trigonometric functions in the collinearity need be evaluated 
only once for each flight line section. It should be emphasized 
that the conclusions drawn and recommendations made are valid only 
for the two flight lines tested. Additional testing of arrays 
generated under varying conditions is necessary before firm conclusions 
may be drawn for digital MSS arrays in general. 

Similar to the collinearity, the piecewise polynomial formulation 
represents another "parametric" method suggested for use in resti- 
tuting MSS digital data arrays. In the next section, test cases 
for this method and results of their application to the data arrays 


chosen will be presented. 
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1+.4.2 Piecewise Polynomials 

To test the suitability of the use of piecewise polynomials 
in analysis and restitution of digital MSS data arrays, the polynomial 
forms of Section 3.2.2 were used. The data arrays were not resampled, 
but were used directly as they came from the digitizing equipment. 

As with the collinearity formulation, each polynomial form chosen 
was tested as the mathematical model, treating the entire flight 
line as one unit, in two approximately equal sections, and in three 
approximately equal sections. The constraints of Equations (3.38) 
were invoked at section interfaces* The arrays were roll stabilized, 
so that a) = 0 was assumed, and the effect of <f> was assumed incorporated 
within the X c variation. 

As previously mentioned in Section 3.2.3, the formulation of 
piecewise polynomials will yield results identical to the collinearity 
equations if the trigonometric series for tan 0 is taken to a sufficient 
number of terms, and if the coefficients of the polynomials are 
estimated using the same least squares models as was used to estimate 
the parameters of the collinearity equations. 

As a test of this supposition, the case was investigated in 
which the orientation elements were assumed to be approximated by 
second order polynomials for the entire length of flight line 218. 

This would represent Case C2, with the flight line as a single unit, 
in the previous section. Polynomials were generated similar to 
those of (3.20) except that tan 0 was represented using additional 
terms in the series expansion beyond the first two represented in 
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(3.20). This formulation, was programmed using the least squares 
model of combined observations and parameters. The a priori 
estimates of variance for the x, y, X, Y observed quantities were 
assigned the same values as were assumed for the collinearity 
formulation. The coefficients of the polynomials were then estimated 
using successive relinearizations in an iterative least squares 
solution. The results obtained were a reference variance of 1.08 
and a positional check variance of 5-04. These values are identical 
to those obtained for the corresponding collinearity analysis, as 
may be seen in Table 4.8. In fact, the results were identical in 
all respects to those from the collinearity analysis. Thus, if 
sufficient terms are taken in expanding tan 0, and if the same rigorous 
least squares procedure is used to estimate coefficients, the results 
of a piecewise polynomial formulation will be identical to those 
using the collinearity equations. 

However, the polynomials in Section 3. 2. 2.1 were derived assuming 
only the first two terms in the series expansion of tan 0. Further, 
it is customary when utilizing polynomials to solve for the coefficients 
by least squares using the method of indirect observations (multiple 
regression). This approximation requires no iterative procedure. 

The estimates of the polynomial coefficients result from a 
single solution of the linearized equations (V + BA = F°) and 

the observations are assumed to have identity weight matrix. These 
approximations were made for' the same test case from flight 218 
and resulted in a positional check variance of 5.83 and an 
a posteriori reference variance from the adjustment of 7 * 78 . 
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The positional check variance for the more rigorous adjustment was 
5.04. The reference variance resulting from the more rigorous 
adjustment was 1 . 08 , which when multiplied by the a priori variance 
estimate of 6.25 assumed for positions in the data array, yields 
a value of 6.75. The differences between these variances from the 
more approximate procedure (5*83 and 7*78) and the rigorous procedure 
( 5 - OU and 6.75) are not statistically significant at an a level of 
0.05* It may reasonably be assumed then that the approximate procedure 
is suitable for use with piecewise polynomials. In the subsequent 
tests, the approximate polynomials of Section 3. 2. 2.1 were used, 
and coefficients were estimated using the least squares model of 
indirect observations. 

The polynomial forms chosen to be tested were the following: 

Case PI: The polynomials for this case were of the form: 


X = A 1 + A 2 x + A^y + Ai + xy + (A^/3c 2 )y 3 + (Al t /3c 2 )xy 3 
Y ■ B-l + B 2 x + B 3 y + xy + (Bg/3c 2 )y + (B^/3c 2 )xy 3 


( 4 . 7 ) 


which represents one of the simplest polynomial forms possible. 

In this case orientation element variations are assumed linear, 
and the first two terms of the trigonometric series expansion are 
used to approximate the scan angle. Topographic effects are neglected 


in this formulation. 
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Case P2: For this case, the polynomials used were: 

X = A 1 +Agx + A 3 y + Aj l ay + (A 3 /3c 2 )y 3 + (Aj | /3c 2 )xy 3 

- Z[(A 3 /Z c )y + (A 3 /3c 2 Z c )y 3 + (A^/Z Q )xy + (Aj + /3c 2 Z c )xy 3 ] 

Y = B 1 + B 2 x + B 3 y + xy + (B 3 /3c 2 )y 3 + (Bj + /3c 2 )xy 3 

- Z[(l/c)y + (l/3c 3 )y 3 ] 

This formulation represents linear modeling of exterior orientation 
variation, but includes the effect of topographic elevations. 

It is equivalent to Case Cl of the collinearity formulation, in that 
the same assumptions concerning exterior orientation variations were 
made. 

Case P3: The polynomials 

X - A^ + A£ x + A^ x 2 + A^ y + A^ xy + Ag x 2 y 

+ (A^/3c 2 )y 3 + (Ac/3c 2 )xy 3 + (A^/3c 2 )x 2 y 3 

* (U.9) 

Y = Bj + B 2 x + B 3 x 2 + Bj^ y + B^ xy + Bg x 2 y 

+ (B^/3c 2 )y 3 + (Brj/3c 2 )xy 3 + (B£/3c 2 )x 2 y 3 
were used. This formulation represents the assumption of second 


order polynomials for the variations of all of the elements of exterior 
orientation, in which terrain elevations are neglected. 
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Case P4: The polynomials used for this case are: 

X = A-^ + x + Ag x 2 + y + xy + Ag x 2 y 

+ (A4/3c 2 )y 3 + (Atj/3c 2 )xy 3 + (Ag/3c 2 )x 2 y 3 

- Z[(A^/Z c )y + (A u /3c 2 Z c )y 3 + (A 5 /Z c )xy 

+ (A 5 /3c 2 Z c )xy 3 + (A 6 /Z c )x 2 y + (A 6 3c 2 Z c )x 2 y 3 ] ^ 

Y = B^ + B 2 x + B3 x 2 + B4 y + xy + B^ x 2 y 

+ (B^/3c 2 )y 3 + (B^/3c 2 )xy 3 + (Bg/3c 2 )x 2 y 3 

- Z[(l/c)y + (l/c 3 )y 3 ] 

This form of the polynomials assumes second order variations of 
orientation elements and attempts to include topographic effects. 

It represents the same assumptions concerning exterior orientation 
behaviour as were made for the collinearity model in Case C2. 

Table 4.13 summarizes the results of the analysis for flight 208. 
Table 4.l4 depicts a tabulation of F statistics for this flight 
line formed from the reference variances from least squares analysis, 
and Table 4.15 is a tabulation of F statistics for this flight based 
upon check variances computed using check points withheld from the 
adjustment . 

It may be seen from these tabulations that for this flight 
over flat terrain, the inclusion of elevation data had virtually 
no effect upon the accuracy obtained, as would be expected. For the 
lower order polynomials of Cane PI and Case P2, statistically significant 
improvement was noted in dividing the strip into 2 sections over 
that of treating the flight line as a unit, and little further 
improvement was noted for these lower order polynomials in going 
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Table 4.13 





Statistics from Piecewise Polynomial Analysis 

of Flight 208 




(60 Check Points) 







Number 

of Sections 




1 

2 

3 

Case 

PI: 

Number of Parameters 

8 

16 

24 



Degrees of Freedom 

TO 

66 

62 



Reference Variance 

7.91 

3.29 

2.87 



Check variance x in rsu 2 (d.o.f. = 59) 

9.10 

3.06 

2.42 



Check variance y in rsu 2 (d.o.f, = 59) 

4.54 

3.96 

3.39 



Positional Check Variance in rsu 2 

6.62 

3.50 

2.89 

Case 

P2 : 

Number of Parameters 

8 

16 

24 



Degrees of Freedom 

70 

66 

62 



Reference Variance 

7.92 

3.29 

2.86 



Check Variance x in rsu 2 (d.o.f. = 59) 

9.10 

3.06 

2.42 



Check Variance y in rsu 2 (d.o.f. = 59) 

4.57 

4.01 

3.41 



Positional Check Variance in rsu 2 

6.64 

3.52 

2.89 

Case 

P3: 

Number of Parameters 

12 

24 

36 



Degrees of Freedom 

66 

58 

50 



Reference Variance 

3.16 

2.41 

1.82 



Check Variance x in rsu 2 (d.o.f. = 59) 

2.49 

2.28 

2.01 



Check Variance y in rsu 2 (d.o.f. =59) 

4.05 

2.84 

1.88 



Positional Check Variance in rsu 2 

3.22 

2.56 

1.95 

Case 

P4: 

Number of Parameters 

12 

24 

36 



Degrees of Freedom 

66 

58 

50 



Reference Variance 

3.16 

2.40 

1.83 



Check Variance x in rsu 2 (d.o.f. = 59) 

2.49 

2.28 

2,01 



Check Variance y in rsu 2 (d.o.f. = 59) 

4.09 

2.86 

1.92 



Positional Check Variance in rsu 2 

3.24 

2.57 

1.96 



Table 4.l4 


Tabulation of F Statistics for Flight 208 
Formed from Reference Variance from Piecewise Polynomial Analysis 



Case PI 
12 3 

Case P2 

Humber of 
12 3 

Case P3 

Sections 

12 3 

Case P4 
12 3 

Case PI 1 

2 
3 

. 

1.00 

2.40* 1.00 
2.76* 1.15 1.00 




n 

Case P2 ° 1 

+» p 
o c - 
4) -3 

CQ 3 

1.00 2.40* 2.76* 
2.40* 1.00 1.15 
2.76* 1.15 1.00 

1.00 

2.41* 1.00 

2.77* .115 1.00 



<M 

0 

Case P3 1 

*4 2 

<u c 

1 3 
a 

2.50* 1.04 1.10 
3.28*1.36 1.19 
4.35* 1.81* 1.57* 

2.51* 1.04 1.10 
3.29* 1.36 1.18 
4.35* 1.81* 1.57* 

1.00 

1.31 1.00 
1.74* 1.32 1.00 


Case P4 1 

2 
3 

2.50* 1.04 1.10 
3.30* 1.37 1.20 
4.32* 1.80* 1.57* 

2.51* 1.04 1.10 
3.30* 1.37 1.19 
k.33* 1.80* 1„56 

1.00 1.31 1.74* 
1.32 1.00 1.32 
1.73* 1.32 1.00 

1.00 

1.32 1.00 
1.73* 1.31 1.00 


* Denotes values which are statistically significant at an a level of .05. 


Table 4.15 


Tabulation of F Statistics for Flight 208 
Formed from Positional Check Variances from Piecewise Polynomial Analysis 

Case PI Case P2 Case P3 Case P4 

Number of Sections 



12 3 

12 3 

12 3 

12 3 

Case PI 1 

2 
3 

1.00 

1.89* 1.00 

2.30* 1.21 1.00 




o 

cn 

o 

ro 

Sections 
oo ro h 

1.00 1.69* 2.30* 

1.88* 1.01 1.22 
2 . 29 * 1.20 1.01 

1.00 

1.88* 1.00 
2.30* 1.21 1.00 



o 

Case P3 1 

h 2 

<j) 

| 3 

2.05* 1.09 1.11 
2 . 60 * 1.37 1-13 
3.41* 1 . 80 * 1.48 

2.05* 1.09 1.11 
2 . 60 * 1.38 1.14 
3.42* 1.81* 1.49 

1.00 

1.26 1.00 
1.65* 1.31 1.00 


Case P4 1 

2 
3 

2.04* 1.08 1.12 
2.59* 1.36 1.12 
3.37* 1.79* 1.47 

2.04* 1.09 1.11 
2.59* 1-37 1.13 
3.37* 1 . 80 * 1.48 

1.01 1.27 1.66* 
1.25 1.00 1.32 
1.64* 1.31 1.00 

1.00 

1.26 1.00 
1.65* 1.31 1.00 


ro 

4T- 

U> 


* Denotes values which are statistically significant at an a level of .05. 
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from 2 to 3 sections. For the higher order polynomials, it was 
in general necessary to segment the flight into 3 sections before 
significant improvements could be discriminated over treatment as 
a single unit. In most cases, significant improvement was noted 
in going from the lower order polynomials of Case PI and Case P2 
to the higher order polynomials of Case P3 and Case P4 for the same 
number of sections . 

Table 1+.16 summarizes the results of the polynomial analysis 
for flight 218. Tables ^.17 and ^.18 depict tabulations of F statistics 
formed from the reference variances and check variances, respectively. 
Although improvements in the reference variances were generally 
noted when elevation information was included for this flight over 
moderate terrain relief, such improvements were not found to be 
statistically significant. Segmenting of the strip had no significant 
effect upon the accuracies obtained for the lower order polynomials. 

The only generally apparent statistically significant improvement 
in the reference variances appears to be in going from the lower 
order polynomial cases to the higher order Case P4 combined with 
segmenting the flight line into 3 sections. No significant differences 
were detected for flight 218 based upon the check variances, as may 
be expected due to the small degrees of freedom associated with 
these (8). Note that for this analysis, as with the collinearity 
model, the assumption of a large number of parameters could be 
dangerous, as indicated by the increase in the check variances. 

Thus, although the points included in the adjustment were fit quite 
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Table 4 . l 6 


Statistics from Piecewise Polynomial Analysis of Flight 218 

(9 Check Points Only) 

Number of Sections 




1 

2 

3 

Case PI: Number of Parameters 


8 

16 

24 

Degrees of Freedom 


38 

34 

30 

Reference Variance 


10.20 

7.94 

6 . 80 

Check Variance x in rsu 2 (d.o.f. = 

8) 

T .20 

6.61 

7.33 

Check Variance y in rsu 2 (d.o.f. = 

8) 

5.67 

4.74 

1.84 

Positional Check Variance in rsu 2 


6.4o 

5.64 

4.13 


Case P2: Number of Parameters 


8 

16 

24 

Degrees of Freedom 


38 

34 

30 

Reference Variance 


9.11 

7.02 

6.33 

Check Variance x in rsu 2 (d.o.f. = 

8 ) 

7.19 

6.58 

7.26 

Check Variance y in rsu 2 (d.o.f. = 

8 ) 

4.96 

4.72 

2.36 

Positional Check Variance in rsu 2 


6.02 

5.61 

4.47 


Case P3: Number of Parameters 


12 

24 

36 

Degrees of Freedom 


34 

26 

18 

Reference Variance 


8.88 

5.52 

4.19 

Check Variance x in rsu 2 (d.o.f. * 

8 ) 

6.67 

14.96 

17.71 

Check Variance y in rsu 2 (d.o.f. - 

8 ) 

4.91 

2.02 

1.32 

Positional Check Variance in rsu 2 


5.76 

6.99 

7.17 


Case P4: Number of Parameters 


12 

24 

36 

Degrees of Freedom 


34 

26 

18 

Reference Variance 


7.78 

5.07 

3.56 

Check Variance x in rsu 2 (d.o.f. - 

8 ) 

6.64 

14.88 

17.69 

Check Variance y in rsu 2 (d.o.f. « 

8 ) 

5.02 

2.13 

1.4o 

Positional Check Variance in rsu 2 


5.81 

7.07 

7.26 



Table 1+.17 


Tabulation of F Statistics for Flight 218 
Formed from Reference Variance from Piecewise Polynomial Analysis 

Case PI Case P2 Case P3 Case Pl+ 


Number of Sections 



12 3 

12 3. 

12 3 

12 3 

Case PI 1 

2 
3 

1.00 

1.28 1.00 

1.50 1.17 1.00 




01 

Case P2 1 

-P 2 
o el 

d> q 

CD ° 

1.12 l.l4 1.34 
1.1+5 1.13 1.03 
1.6l 1.25 1.07 

1.00 

1.30 1.00 
1.1+1+ 1.11 1.00 



0 

Case F3 1 

u 2 
<u c - 

1 3 

1.15 1-12 1.30 
1.85* 1.1+1+ 1.23 
2.1+3* 1.89 1.62 

1.03 1.26 1.1+0 
1.65 1.27 1.15 
2.17* 1.68 1.51 

1.00 

1.61 1.00 
2.12* 1.32 1.00 


Case PU 1 

2 
3 

1.31 1.02 1.1k 
2.01* 1.56 1.34 
2.86* 2.23* 1.91 

1.17 1.11 1.23 
1.80 1.38 1.25 
2.56* 1.97 1.78 

1.11+ 1.1+1 1.86 
1.75 1.09 1.21 
2.1+9* 1.55 1.17 

■ 

1.00 

1.53 1.00 
2.18* 1.1+2 1.00 


ro 

C\ 


* Denotes values which are statistically significant at an a level of .05. 


Table 4.18 


Tabulation of F Statistics for Flight 218 
Formed from Positional Check Variances from Piecewise Polynomial Analysis 

Case PI Case P2 Case P3 Case P4 

Number of Sections 


' 

12 3 

12 3 

12 3 

12 3 

Case PI 1 

2 
3 

1.00 

1.13 1.00 
1.55 1.3T 1.00 


t 


« 

Case P2 .2 1 

-p p 
o eL 
o 

to J 

1.07 1.06 1.46 
l.lU 1.01 1.36 

1.43 1.26 1.09 

1.00 

1.07 1.00 
1.35 1.25 1.00 



0 

Case P3 1 

2 

a; <L 

1 3 

1.11 1.02 1.39 
1.09 1-24 1.69 
1.12 1.27 1.74 

1.04 1.03 1.29 
1.16 1.25 1.56 

1.19 1.28 1.60 

1.00 

1.21 1.00 

1.25 1.02 1.00 


Case P4 1 

2 
3 

1.10 1.03 1.4l 

1.10 1.25 1.71 

1.13 1.29 1.76 

1.03 1.04 1.30 
1.17 1.26 1.58 
1.21 1.29 1.62 

1.01 1.21 1.24 
1.23 1.01 1.02 
1.26 1.04 1.00 

1.00 

1.22 1.00 

1.25 1.02 1.00 


ro 
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closely as shown by the lower reference variances, the true test 
of the fit, as given by the check variances computed from points 
withheld from the adjustment, deteriorated. 

The value in including Z coordinate elevations is doubtful 
for this test array over moderate terrain relief. For this test, no sig- 
nificant improvement was gained in doing so, indicating that terrain 
relief effects are relatively small for this array, and such effects 
are at least partially compensated for in estimating the parameters 
of exterior orientation. However, a test with many more check points 
should be carried out before a more concrete assertion can be made 
regarding the inclusion of point elevations in the polynomials . 

4.4.3 Nonparametric Methods 

The four nonparametric methods of MSS digital data array 
analysis and restitution which were attempted for testing are: 

1. Arithmetic mean 

2 . Moving average 

3. Meshwise linear transformation 

4. Linear least squares filtering 

The first three of the above represent algorithms which may be 
applied to MSS digital data arrays in a fairly straightforward 
manner. The linear least squares filtering is highly dependent 
upon the validity of assumption of some statistical properties 
concerning the data arrays. Considerable difficulty was experienced 
in applying this method to the data arrays available for testing. 

In implementing the arithmetic mean algorithm, the data arrays 
were first preprocessed using the resampling algorithm of 
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Equation (3-13), followed by approximate scaling in the y direction, 
in order to approximately equalize scales in the x and y directions. 
A weighting function was then chosen using only separation distance 
between points as the independent variable. For the flight lines 
tested, the weighting function was assumed to be of the form 


1 

w i = W (k.H) 

For the test of the moving average algorithm, no preprocessing 
of the data was done. The weighting function of (U.ll) was also 
used for this method, and the non-linear form of Equation (3.1+5) 
was utilized. 

The meshwise linear transformation was also tested on both 
flight lines. No preprocessing of the data is required for this 
method, and no weighting function need be defined. The control 
points were connected into a triangular mesh, and each check point 
was transformed according to the unique parameters of an affine 
transformation defined by the coordinates of the three control points 
forming the triangle within which the check point lies. 

Attempts were also made to analyze the applicability of the 
linear least squares filtering method as a restitution technique. 

The results of these attempts, however, were disappointing. Suc- 
cessful application of the method depends upon the data having 
specific statistical properties. The property of ergodicity must 
be assumed, which in turn depends upon stationarity of the data. In 
addition, isotropy of the two dimensional data field must be assumed [99]- 
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The resampling algorithm of Equation (3.13), followed by ap- 
proximate y direction scaling was used to preprocess the flight 
lines in an attempt to satisfy the isotropic requirement. Sample 
covariance points were then computed from the data sets using the 
method presented in [993* Results of these steps for flight 208 
are represented by Figure 4.11. From these plots, it is apparent 
that the data set cannot be considered to satisfy the statistical 
requisites over the entire length of the flight line, since no logical 
fit to these sample covariance points may be obtained using an allowable 
autocovariance functional form. Yaglom [1093 and Bendat and Piersall 
[114] enumerate the characteristics of an autocovariance function 
as follows : 

C(0) > 0 

C(-x) = C(t) 

C(t) - C(0) 




XY Covariance (C xy ) Y Covariance (C y ) X Covariance (C^) 


25] 



Figure i+.ll. Sample Covariance Points for Flight 208 
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forms, as suggested by Kraus and Mikhail [99], to sample covariance 
points. The Gaussian type covariance function may be written in 
the form: 

C(d) = C(0) e“ kd2 (4.12) 


in which the C(0) and k parameters are estimated by a least squares 
fit of the function to the sample covariance points . 

After these covariance functions were determined, they were 
used to calculate the elements within the C_ matrix of Equation (3.49), 
as well as the c_ matrix. Each element of the sub-matrices C x , Cy, 

Cxy, _Cy X in Equation (3.49) are computed from the appropriate covariance 
function. For the present investigation, two cases were considered. 

In the first, three separate covariance functions were determined 
for calculation of the terms in the C x , Cy and C^y sub-matrices 
by a least squares fitting to the sample covariance points. For 


this case it was assumed C^y = Cy X . In the second case. 


it was 


assumed that C^y = _C yx = 0. This represents the case for which 

independent interpolation and filtering are done in the x and y 


directions . 


If interpolation only is to be done, the diagonal elements 
of the £ matrix will contain values C(0) from the appropriate auto- 
covariance functions. However, if estimates of random components 
r at the reference points are to be filtered, in addition' to interpo' 
lation of the s components, the diagonal terms of the jC matrix are 
replaced by variance terms V x , V y computed from the data set. 
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For the x variance terms, for example, the variance is computed 

by 


V x = Z 4 2 xi 
i-1 

in which £ . = 3L - for n control points. It was this technique 

which was attempted in this investigation. Hence, simultaneous 
interpolation and filtering was attempted. 

For the data sets available for testing, however, the results 
obtained using this linear least squares procedure were highly 
erratic. Combinations of covariance functions which appeared to 
yield acceptable results on one flight line yielded very large 
positional check variances on the other. A preliminary supposition 
is that the data sets, even after the preprocessing steps described 
above, do not adhere to the requisite statistical properties necessary 
for the successful application of the linear least squares algorithm. 

Because of this unpredictability of results from the linear 
least squares procedure, this technique is not included in the 
tabulated results to be presented. In Table 4.19 the resulting 
statistics obtained from the processing of the data arrays for flight 208 
by nonparametric interpolative methods are shown. The statistics 
shown are based upon check variances formed by applying the non- 
parametric methods to the check points as based upon the control 
points available within the strip. Sixty such check points were 
available for flight 208. In Table 4.20 a similar tabulation is 
shown for flight 218. Only nine check points were available for this 
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Table 4.19 


Statistics from Analysis of Flight 208 
by Nonparametric Interpolative Methods 
(60 Check Points) 


Arithmetic Mean 


Degrees of Freedom (x and y) 
Check Variance, x, in rsu 2 
Check Variance, y, in rsu 2 
Positional Check Variance 


59 

2. kb 
1.51 
1.86 


Moving Average 

Degrees of Freedom (x and y) 
Check Variance, x, in rsu 2 
Check Variance, y, in rsu 2 
Positional Check Variance 

Meshwise Linear Transformation 


Degrees of Freedom (x and y) 
Check Variance, x, in rsu 2 
Check Variance, y, in rsu 2 
Positional Check Variance 


59 

1.75 

4.16 

2.79 


59 

1.82 

5.09 

3.26 
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Table U.20 

Statistics from Analysis of Flight 218 
by Nonparametric Interpolative Methods 
(9 Check Points) 


Arithmetic Mean 


Degrees of Freedom (x and y) 
Check Variance, x, in rsu 2 
Check Variance, y, in rsu 2 
Positional Check Variance 


Moving Average 

Degrees of Freedom (x and y) 
Check Variance, x, in rsu 2 
Check Variance, y, in rsu 2 
Positional Check Variance 


Meshwise Linear Transformation 


Degrees of Freedom (x and y) 
Check Variance, x, in rsu 2 
Check Variance, y, in rsu 2 
Positional Check Variance 


8 

10.69 

2.32 

5.76 


8 

7 . 1*9 

3.82 

5.52 


8 

6.25 

5.86 

6.05 
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flight line, and the resulting degrees of freedom (8) for the x and 
y coordinate values, respectively, is quite small. Table 4.21 is 
a tabulation of the F statistics computed from the positional check 
variances for flight 208. Of interest here is the fact that the 
simplest interpolative method, that of the arithmetic mean, shows 
a significant reduction of variance over the meshwise linear transfor- 
mation method. No other statistically significant differences were 
noted for this flight line. Table 4.22 shows a similar tabulation 
for flight 218. As might be expected, no significant differences 
between the interpolative methods were detected for the small degrees 
of freedom available. The inference to be drawn from these tables 
is that if restitution is to be done for the arrays tested using 
a nonparametric method, the choice of the simple arithmetic mean 
represents a logical and computationally efficient choice , well 
suited to automated restitution as a portion of automated interpretation 
General conclusions should not be drawn, however, based upon such 
limited tests. Further investigations utilizing data arrays generated 
under a variety of conditions must be performed before a definitive 
statement may be made. 

4.4.4 Summary of Results 

The results of the tests performed on the data arrays chosen 
indicate that the nonparametric methods, particularly the arithmetic 
mean algorithm, represent perhaps the optimum method for restitution 
of MSS digital data arrays. For flight 208, each of the check var- 
iances from the collinearity and piecewise polynomial analysis were 
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Table 4.21 

Tabulation of F Statistics for Flight 208 
Formed from Positional Check Variances from Nonparametric Methods 



Arithmetic 

Mean 

Moving 

Average 

Meshwise 

Linear 

Trans formation 

Arithmetic 

Mean 

1.00 



Moving 

Average 

1.50 

1.00 


Meshwise 

Linear 

Trans formation 

1.75* 

1.17 

1.00 


* Denotes values which are statistically significant at an a level 
of .05. 


Table 4.22 

Tabulation of F Statistics for Flight 218 
Formed from Positional Check Variances from Nonparametric Methods 



Arithmetic 

Moving 

Meshwise 


Mean 

Average 

Linear 

Tran s formation 

Arithmetic 

Mean 

1.00 



Moving 

Average 

1.04 

1.00 


Meshwise 

Linear 

Trans formation 

1.05 

1.10 

1.00 
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compared using the F statistic to that obtained by the arithmetic 
mean algorithm. It was found that, for the collinearity analysis, 
only the higher order polynomial assumptions, coupled with dividing 
the flight line into sections could compare favorably. In all cases 
of linear orientation variation assumption the check variances from 
the collinearity analysis were significantly higher than that obtained 
from the arithmetic mean analysis, at an a level of .05. Similar 
results were obtained when a comparison of the results of the poly- 
nomial analysis wi'‘:h that of the nonparametric methods was done. 

For flight 208, only the higher order polynomials coupled with 
dividing the flight line into three sections could compare favorably 
with any of the nonparametric methods . The nonparametric methods 
in general produced statistically better results than those from 
the lower order polynomials. 

For flight 218, a comparison of the check variances obtained 
from the collinearity analysis and those from the nonparametric 
methods was also done. No statistically significant differences 
were found. A comparison of check variances from the polynomial 
analysis with those from the nonparametric methods also revealed 
no significant differences. These results are noteworthy, because 
in using the nonparametric methods no attempt was made to account 
for displacements due to topography. This would indicate that for 
arrays obtained by scanning areas of moderate terrain relief (e.g. 
flight 218) that errors due to relief displacement are small compared 
to those due to other sources (orientation, scan time, etc.), or 
that the nonparametric algorithms at least partially account for 
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the topographic effects or a combination of these factors. There- 
fore the nonparametric methods (particularly the relatively simple 
arithmetic mean algorithm) were shown to represent an effective 
method for use in the restitution of such arrays, at least 
for the arrays tested, although no general conclusions should be 
drawn from such limited tests. 

4.5 Final Transformation of Flight Lines 
A final transformation of the test flight lines was attempted 
utilizing the collinearity method of restitution. The method used 
was as described in Section 3.5. After estimating the parameters 
of the model, the restituted Xj , Yj coordinate values were assigned 
as the independent variables in *,he transformation equations. An 
iterative procedure was utilized in which Zj was first estimated, 
the collinearity equation for the coordinates in the direction of 
flight was then solved explicitly for Xj, the uncorrected array 
position, and this value was utilized in the second collinearity 
equation to solve for yj, the along scan array position in the 
uncorrected array. These values were used to locate an updated 
value for the Zj term, and the process was repeated until no change 
in the resulting array position was found. The spectral values in 
the Xj , y^ uncorrected array position were then assigned to the 
Xj, Y^ restituted array position after iterative transformation. 

This procedure was repeated for each data array element. The program 
utilized was a modification of Anuta's OVERLA program [110], which 
was originally written to superimpose different arrays of data con- 
taining essentially the same ground area. 




260 


Digital displays of the restituted arrays were then generated 
to assess qualitatively the results. In Figure 4.12 is shown the 
digital display of channel 6 for flight line 208 after restitution, 
which may he compared directly to Figure 4.2, the digital display 
of the unprocessed array. A careful study of these images reveals 
that improvements in the geometric characteristics of the array 
have teen obtained. The narrowing of the flight line results from 
a correction of the scan angle effects and the differential scale 
in the x and y directions. The triangular area near the top of the 
flight line (arrow A, Figure 4.12) has been visibly restored to a 
configuration more closely matching the shape it should have as 
illustrated in the map planimetry of Figure 4.5. The alignment of 
the interstate highway cutting diagonally across the bottom half 
of the flight line (arrow B, Figure 4.12) has been visibly straightened, 
as has the alignment of the central roadway running down the flight 
line. Cross roads appear more nearly perpendicular to this central 
roadway, indicating that the scan time effects and the effects of 
aircraft yaw were compensated for. 

The restituted array from channel 6 for flight line .218 is 
shown in Figure 4.13. The river bend (arrow A, Figure 4.13) in 
this display is seen to reflect more closely the map planimetry 
illustrated in Figure 4.6 than in the unprocessed array of Figure 4.4. 
The inclination of the roadway (arrow B, Figure 4.13) in this area 
also has been restituted into an alignment more closely matching 
that shown in Figure 4.6. 
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These displays illustrate the feasibility of algorithms to 
restitute the originally deformed data arrays into array systems 
displayable in an orthographic projection system. These restituted 
data arrays may be used to extract metric as well as interpretive 
information, and may serve as a basic data source for the existing 
automated interpretation algorithms , reflecting more closely the 
geometry of the area sensed as it actually appears in the ground 


system. 



5. CONCLUSIONS AND RECOMMENDATIONS 


The scope of the investigation has encompassed the analysis 
of the geometric aspects of recording multispectral scanner (MSS) 
data in a digital manner. An introduction to the digital techniques 
of utilizing MSS data has been presented, and the geometric distortions 
present in such arrays were described. The basic transformation 
equations for recording MSS data were given and specific expressions 
were derived from these for analysis and restitution of the data 
arrays. Nonparametric techniques of restitution were also utilized 
as an alternative to the parametric methods based upon the trans- 
formation relationships. Testing of both methods of analysis and 
restitution on actual digital data arrays was attempted. The 
numerical results obtained were for the case of singly scanned flight 
lines only, since no overlapping flight lines were available for 


testing in this investigation,. 

\ 

X 

5.1 Con cl 


us ions 


The investigation conducted supports the following conclusions: 
1.) Significant geometric distortions are present in unprocessed 
fCS digital data arrays, and must be considered in order to reliably 
obtain reliable metric information for surveys^ of resources at 


\ 

\ 
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the Earth’s surface. The major causes of these geometric pertur- 
bations present within the arrays are: 

a. ) changing ground resolution size at different scan angles 

and image position displacement due to panoramic recording 

b. ) terrain relief displacements due to topographic relief 

differences within the area scanned 

c. ) perturbations caused by the continuous changes in sensor 

exterior orientation elements during recording of the data 

d. ) effective image displacements due to forward motion of the 

sensor during the short time interval required to record 

a single scan line. 

Two types of investigations were carried out for data arrays containing 
such geometric displacements. The problem of analysis was considered, 
in which it is desired to investigate the causes of the geometric 
displacements, and determine the resulting magnitudes and directions 
of the errors within the data arrays due to each cause. Subsequently, 
the problem of restitution of the data arrays was considered, in 
which it is desired to correct the geometric errors within the arrays 
using any convenient method, without necessarily isolating and 
investigating the individual causes of the errors. 

2.) The MSS collinearity equations derived from the geometric 
portion of the more general MSS transformation provide a suitable 
mathematical model for both analysis and restitution of MSS data. 

These expressions represent the actual geometric relationships 
between ground points and image array positions at the moment of 
recording. Therefore they represent the most obvious mathematical 
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model, and serve as the basis of other approximate methods of analysis 
and restitution. 

3. ) The collinearity equations may be used with ground control 
points to perform space resection in order to obtain estimates for 
the elements of exterior orientation of the sensor during recording 
of the data. However there are some severe limitations on such 
resections. The most serious is that, since the multispectral scanner 
records continuously, each scan line has different exterior orientation 
elements than any other scan line. To solve for a separate set 

of exterior orientation elements for every scan line is impossible. 
Therefore, some assumption must be made concerning the variation 
of each exterior orientation element as the sensor progresses down 
the flight line. If a functional form is assumed for each of these 
orientation elements, the parameters of these functions will be de- 
termined by a space resection. If the actual stochastic variations 
of the exterior orientation elements are not reasonably approximated 
by the functional forms assumed, then the collinearity formulation 
and methods based upon the collinearity formulation lose their ad- 
vantage, and other methods of analysis may compete favorably with 
tha collinearity formulation as a mathematical model. It was found 
in this investigation that the assumption of polynomial functional 
forms for the orientation elements provided adequate accuracies 
and reasonable computational efficiency foi 1 the data arrays tested. 

4. ) If space resection is performed utilizing ground control 
points based upon the MSS collinearity equations for a singly scanned 
flight line, not all of the elements of exterior orientation for the 
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sensor can be recovered. A constant pitch ($) of the sensor is 
not distinguishable from a constant X translation of the sensor 
perspective center (X c ). Likewise, a constant rate of change of 
pitch is not distinguishable from a constant aircraft velocity. 

Similar dependencies are found if higher order terms are attempted 
to simultaneously model the pitch and X translation of the sensor. 
Therefore not all of the elements of exterior orientation are re- 
coverable . 

5. ) Significant improvements in the results of restitution 
using collinearity equations may be realized by breaking the flight 
line into sections. For the flight lines tested, sections of 800 
lines or less were used. Control points within each section can 
then be used to define independent parameter sets for that section 

by space resection. Discontinuities in the restituted array positions 
at section boundaries may be avoided by utilizing parameter constraints 
during. the least squares estimation of the exterior orientation 
parameters. This problem of discontinuities at section boundaries 
will be discussed further in a later conclusion. 

6. ) If restitution of MSS digital data arrays is to be attempted 
using the collinearity equations, then terrain height information 
may be introduced into the uncorrected data arrays in some manner, 

in order that relief displacements may be removed during the restitution 
procedure. One possible source of this information (used for this 
investigation) is the U.S. Geological Survey 7 1/2 minute quadrangle 
topographic map sheets. If these or other topographic map data 
are used as a source of elevation information, the data contained 
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on the map must "be transformed into a digital data format in some 
manner. This digitized data must be superimposed upon the uncorrected 
digital MSS arrays such that an elevation is assigned to each element 
within the ftES digital data array with sufficient accuracy to allow 
restitution of the arrays to a precision within the size of the 
resolution of the data. 

7. ) Piecewise polynomial forms for restitution of MSS digital 
data arrays may be derived from the collinearity equations. These 
polynomials may be formulated to account for various assumptions 
concerning exterior orientation element variations, and for relief 
displacements within the arrays. These polynomial formulations 
require that a series expansion be used to approximate the trigono- 
metric functions resulting from scan angle variations. If sufficient 
terms are taken in the series approximations, and if the proper 
rigorous least squares procedure is used, the results obtained 

are identical to the results using the collinearity formulations. 

If only a few terms are carried for the series expansion, and a 
simpler least squares procedure is used, the results from these 
polynomial models are in general, inferior to the results from col- 
linearity analysis. The difference, however, was found not to be 
statistically significant. 

8. ) Both the collinearity and piecewise polynomial formulations 
were found to be relatively insensitive to moderate errors in assignment 
of element elevations to the uncorrected data arrays. The derivation 

of an allowable height assignment error formula revealed that an 
element elevation error of 1.2# of the flying height above terrain 
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would result in an error in position of less than the ground reso- 
lution size of an element for arrays digitized at 6 mi Hi radians , 
at the extreme case of a scan angle of U 5 0 from the vertical. In 
conducting tests on real data having flat and moderate terrain relief, 
it was found that the resulting positional error was affected to 
a greater degree by the assumption of functional forms for exterior 
orientation elements than by the inclusion of element elevation 
information within the expressions utilized. 

9. ) The assignment of elevations to element positions within 
the uncorrected arrays to the required accuracy was found to be 
relatively insensitive to the approximate planimetric transformation 
used to transform the original array positions to the coordinate 
system based upon digitized control points. An affine transformation, 
applied to sections of the flight line was found adequate for this 
purpose. Tests of real data arrays revealed that assignment of 
element elevations utilizing the affine transformation resulted 

in element elevation errors less than the theoretical allowable 
height assignment error. 

10. ) The assignment of element elevations by digitally searching 
for contour-scan line intersections followed by linear interpolation 
for elements between these intersection points was found to result 

in acceptable values for element elevations for the purposes of 
restitution. A more general interpolative technique (arithmetic 
mean) was attempted and found to result in much greater computation 
time with no significant improvement in accuracy. 
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11. ) The use of MSS data in digital array form allows for two 
alternative transformation procedures for restitution. In the first 
instance, the more conventional procedure of utilizing image coordinates, 
or, for MSS digital data, array positions as independent variables 

is pursued. These are utilized in the restitution functions to 
generate restituted arrays having improved geometrical properties. 

This method has the advantage that the projective relationships 
may be used directly and requires no iterative steps in restitution. 
However, the method requires a large output buffer, interpolation 
of spectral density data, and may result in multiple assignments 
of spectral values to a single element location. 

The second possible method utilizes the output array positions 
as independent variables and involves a solution of the restitution 
equations for the corresponding uncorrected data array positions. 

An input buffer is required in this case, and an iterative procedure 
is necessary, since elevation data are associated with the uncorrected 
array positions. 

12. ) If discontinuities are to be eliminated at section boundaries, 
constraints may be enforced to accomplish this. If the flight line 

is divided along lines of constant x image coordinate, the general 
constraint to be enforced is that the restituted X, Y position computed 
for every element along this boundary must be the same when computed 
using the parameters for the section on either side of the boundary. 

This results in six constraint equations whether utilizing the 
collinearity or piecewise polynomial formulations. 
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13. ) The nonparametric methods introduced appear to be quite 
promising as an alternative to restitution of MSS digital data arrays 
using collinearity equations or piecewise polynomials. These non- 
parametric interpolative methods have some advantages over the 
more conventional methods of restitution which make them appear to 
be well suited for application to digital MSS data arrays. 

a. ) The accuracy of restitution using these nonparametric methods, 

as indicated by tests on real data arrays having flat and 
moderate terrain relief, was comparable to the restitution 
accuracy achieved by the more conventional methods utilizing 
collinearity equations or piecewise polynomials. 

b. ) The problem of assignment of element elevations may be 

avoided. The treatment of' the geometric displacements 
as a realization of a two dimensional stochastic field 
from which trend has been reasonably removed, appears to 
have at least partially corrected for relief displacements 
within the interpolative algorithms themselves, without 
the need for assignment of element elevations for every 
data array position. 

c. ) Some of the nonparametric algorithms, notably the arithmetic 

mean algorithm, are quite simple and relatively fast com- 
putationally, resulting in efficiency in transforming 
large numbers of array point positions. 

d. ) There is no problem of breaking the flight line into sections, 

and hence no problem of removing discontinuities at section 
boundaries. The nonparametric algorithms are continuous 
throughout the flight line. 
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The primary disadvantage of the nonparametric interpolative 
methods is that some of them depend upon certain assumptions 
(isotropy, stationarity , ergodicity) of the stochastic field. 

MBS digital data arrays must be preprocessed in some manner to 
approximate isotropy. Serious doubts were raised during the inves- 
tigation as to whether the digital MSS arrays obtained from aircraft, 
may be considered stationary, even over limited areas. For most 
of these nonparametric methods, however, the results of restitution 
were found to be acceptable even though the assumptions were not 
fully realized. 

l4.) For the arithmetic mean algorithm, it was found that a 
simple resampling algorithm, followed by approximate scaling in the 
direction normal to flight could be used to acceptably approximate 
isotropy. 

15 •) Of the nonparametric algorithms, the arithmetic mean al- 
gorithm appeared to be most promising for the data arrays tested. 

The method is simple, fast, and the resulting accuracy of restitution 
is comparable or superior to any of the other nonparametric methods 
or the parametric methods ( collinearity equations and piecewise 
polynomials). The lone disadvantage of the method is that it requires 
pre-processing of each array position as discussed in conclusion l4. 

l6.) The moving average interpolative algorithm appears to be 
able to compete with that of the arithmetic mean. The ability of 
this algorithm to include non linear effects allows for the omission 
of pre-processing of the data arrays. However, for each element 
calculated during restitution using this method, a unique parameter 
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set must "be evaluated by least squares techniques, resulting in an 
algorithm which is slower than that of the arithmetic mean, with 
no statistically significant improvement in accuracy of the restitution 
for the arrays tested. 

17. ) The meshwise linear transformation does not appear promising 
as a nonparametric algorithm for restitution of MSS digital data 
arrays. The method is predicated upon the fact that each element 

to be transformed lies within a triangle whose vertices represent 
control points which are identifiable both upon a map sheet and 
within the uncorrected data arrays. Such a condition cannot be 
fulfilled for every data array element within a flight line. A 
further disadvantage of the method is that a large amount of human 
intervention is required, in order to define a unique mesh of con- 
tiguous triangles formed from the control points. Therefore the 
method does not appear well suited to an automated technique combining 
the interpretive and geometric analysis of the data arrays. 

18. ) The results obtained in attempting to utilize the method 
of linear least squares interpolation and filtering were highly 
erratic. The method is highly dependent upon the fulfillment of the 
assumptions concerning stationarity , ergodicity, and isotropy of 
the stochastic field. The preliminary supposition as of this time 
is that one or more of these assumptions were not fulfilled for the 


data arrays tested. 
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5.2 Recommendations 

The following suggestions are offered for further research 
and development efforts in the field of MSS digital data array 
utilization bases upon the experience of this investigation. 

1. ) An attempt should be made to implement a unified approach 
to amalgamate the classically separated activities of interpretation 
and mensuration. The digital method of analysis of MSS data arrays 
appears ideally suited to this purpose. The inclusion of the analysis 
and restitution methods presented in this thesis with existing 
automated interpretation techniques would result in the ability to 
extract quantitative information in addition to the qualitative 
information presently being generated, thus increasing the information 
obtainable from the MSS data concerning Earth resources. 

2. ) An investigation should be carried out concerning the 
applicability of functional forms other than polynomials to approximate 
the stochastic variation of exterior orientation parameters within 
sections of the flight line. Functional forms such as harmonics, 
although difficult to deal with, may result in greater accuracy 

of restitution. 

3. ) An investigation of the actual variation of exterior orientation 
elements with time should be carried out. Ideally, this would be 
accomplished by fully monitoring all sensor exterior orientation 
elements within an accuracy obtainable by resection. However, 

this would be quite expensive and is unlikely. Alternatively, the 
monitoring or stabilization of one additional exterior orientation 
element (in addition to the roll, w) such as the pitch ($) term could 



275 


be attempted. If this monitoring were combined with a flight over 
an area having very dense control, such that a resection for the 
remaining parameters could be carried out at very frequent intervals 
down the flight line, then a reliable estimate of the actual exterior 
orientation elements could be obtained. It would then be possible 
to state with more confidence the actual errors resulting from 
the omission of one of the exterior orientation elements within a 
section. 

U.) It would be most useful to generate some experimental data 
arrays which include an accurate recording of time. In this way, 
the functions assumed for variation of exterior orientation elements 
could be formulated directly as time functions , and the suitability 
of utilizing the image x coordinate as a measure of time could be 
tested. 

5. ) Analysis and restitution methods should be tested on simulated 
data, having known stochastic characteristics. This would allow 

for controlled experimentation and reveal insights into the behaviour 
of real data arrays. 

6. ) All of the methods of restitution presented in this investi- 
gation should be tested on data arrays generated by flights over 
ground areas having more extreme relief differences. Of particular 
interest would be the applicability of the nonparametric algorithms 
to such data, since the algorithms presented do not consider terrain 
elevations. 


\ 

\ 


7.) An investigation of the use of generalized nonparametric algo- 
rithms should be pursued, in which Z coordinate information is included. 
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8. ) An investigation should "be carried out concerning the 
effect upon accuracy of the control point distribution within a 
flight line. This should include studies of how the number and 
configuration of control points available affects accuracy. Results 
may then indicate the minimum number of control points which must by 
available in order to produce restituted data arrays of a certain 
accuracy requirement, and indicate the optimum configuration of such 
control points within the flight line area. 

9. ) Attempts should be made to utilize sources of planimetric 
control points other than maps. One possibility in this area is 
obtaining metric photography of the area as a control source. 

Aerial triangulation utilizing such photography could then be performed 
to define X-Y-Z ground control point coordinates. If such photography 
were available and was digitized, investigations to further improve 
accuracies could be carried out. Digital correlation computation 
could be investigated as a method of reducing the a priori variance 
of assigning array positions to control points . Digital correlators 
using digitized metric photography and the MSS data arrays should 
be developed in order to assign these control point array positions 
with greater reliability than is possible using subjective human 
decisions . 

Another method of control which should be considered is direct 
targeting. In this manner, ground control point coordinates could 
be determined by direct measurement. This would call for an investi- 
gation of suitable target forms for multispectral scanners. Special 
targets would have to be designed for multispectral scanners, which 
would be discemable in several spectral bands. 
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10. ) A further investigation should be conducted into methods 
of assigning elevations to data array positions. The nonparametric 
algorithms and mathematical modeling methods should be emphasized 
in this investigation, in order to see if better accuracies can 

be obtained, particularly for ground areas with extreme relief 
differences . 

11. ) Further investigations should be conducted concerning 

the applicability of the linear least squares interpolation technique 
for restitution of MSS digital data arrays . These investigations 
should include statistical testing of MSS digital data to ascertain 
whether the properties of stationarity , isotropy, and ergodicity 

K 

may properly be assumed for such data. If so, further investigations 
may be conducted concerning pre-processing algorithms and covariance 
functions which may be useful in applying the linear least squares 
method to MSS digital data arrays . 

12. ) A systems analysis should be attempted to try to determine 
the most effective restitution method. In this investigation, 
accuracy 0.7 restitution has been the primary criteria for comparison. 
This should be supplemented by an economic analysis to determine 

the most efficient of the methods in terms of the total computing 
time involved as well as the human intervention necessary. 

13. ) Algorithms should be developed to determine quantitative 
information during the process of automated interpretation and 
restitution. Of particular interest would be algorithms to determine 
areas on the Earth's surface, iti order that a processing system 

would result in both classifications and estimations of areas containing 


classes of interest. 
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lU.) Investigations should he carried out on the usefulness, 

\ 

of multi-scanned data arrays, particularly sidelapping flight lines. 

- \ 

It is necessary that some data arrays he generated utilizing data 

\ \ 

from sidelapping flight lines, in order to assess the usefulness 
of such data. If this data were to become available, it would he 
possible to perform tests to determine whether this method could 
compete favorably with map contour digitization as a source of element 
elevation information, in terms of accuracy and cost. 

As with most research efforts, the investigations carried out 
point- to new problems and possibilities associated with the utiliza- 
tion of MSS data arrays. There is certainly no lack of possible 
research topics to be considered within this relatively young technology. 
As further research is carried forth, new paths of inquiry will 
doubtless become apparent. 
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APPENDIX A 

DIFFERENTIAL FORMULAS FROM LINEARIZED COLLINEARITY EQUATIONS 

If the collinearity equations of (2.19) are linearized using 
the procedures of Section 2.4, the resulting linearized form of 
the equations is 

• •••••* 

F = F'° + B! D 5x + A' 6 y + B' A + B* C A (A.l) 

— — — — — o — — 

in which the B^, A'^, B T , £, D Jacobian matrices and F'°, 4 , 4 
vectors are as defined in Section 2.4. These equations may he 
written in the form 

• •• • • 

F = F'° + B» D fix + A’ 6 y - I ["fix"! + B* 6 Z + B* C 4 (A.2) 

L«J 3 

in which 



and is obtained by partitioning the B' matrix ( B_ f = [-1 

This linearized form of the collinearity equations is useful, 
as it may be solved directly for (<5X, 6Y) object space differential 
changes, which results in the form 



289 


f 6x1 * F'° + B! D 6x + A* 6y + B' 6Z + B* C A (A.4) 

[jSyJ “ _1 ~ 3 


which may be used to investigate object space pianimetric coordinate 
displacements resulting from small changes in the parameters of 
exterior orientation. 

As an example, consider the case of a single scan line for 
which the elements of exterior orientation may be assumed constant, 
due to the short time interval needed to record a single line. 

For this case 


X = X° + 6X 
Y = Y° + 6Y 


(A.5) 


The term h = Z Q - Z is introduced, and the initial approximations 
of (X, Y) may be defined as (using 0 = y/c) 




+ h tan 0 


0 

2,1 


(A.6) 


and the initial approximations of the angular orientation elements 
( a), k ) may be taken as zero (M * I_) . The F’° vector becomes 



x c - x c 

lY c - Y c - h tan0° + h tan0°_ 


0 

2,1 
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Evaluating the Jacobian matrices of Equation (A. 4), 


bi 

2 % 


A'] 

2,1 



1 

0 

0 

0 

-h -h 


_0 

1 

tan0 

h(l + tan 2 0) 

0 ( 


- 

0 



0 

= 




b’ 3 - 



h 

sec 2 0 


2,1 

tan0_ 


_c 






D = 0 

6,1 6,1 


C = I 

6,6 6,6 


Substituting these results into Equation (A. 4) results in 
6X = 6X - h6 4> ~ h tan06« 

(A. 7) 

SY = h/c sec 2 0 6y - tan0 6Z + 6 Y c + tan0 6Z C + h(l + tan 2 0)6w 

If it is assumed that terrain elevation (z) is constant in order 
to accomplish the transformation of a two dimensional vector then 
Z “ 0. If it is desired to investigate the effects of changes 
of exterior orientation elements only, then the image y coordinate 
is assumed to have no error, so that 6y = 0. Making the substitution 
sec 2 0 = 1 + tan 2 0, the resulting equations become 

SX * 6X - h( 6 4> ) - h tan0( 6k) 

(A.8) 

5Y ~ <SY C + tan0(6Z c ) + h sec 2 0(6w) 

If these results and those of Equations (A. 6) are substituted directly 
into Equations (A. 5), the resulting differential form of these 
equations becomes 
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X - X„ + 6X„ - h6<J> - h tan© 6 k 

C ° (A. 9) 

Y - Y„ + h tan© + 5Y + tane SZ + h sec 2 6u} 
c c c 

Equations (A. 8) may be seen to be identical in form to those 
of Equations (2.1 6 ), demonstrating the validity of the linearized 
collinearity equations as a general form from which specific cases 
may be examined. 
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APPENDIX B 

DERIVATION OF RESAMPLING ALGORITHM TO INCLUDE TOPOGRAPHIC EFFECTS [ 96 ] 

Referring to Figure 3.^, the scan line length is given by 

W = hj_ tan(ny) + h K tan(Ny - ny) (B.l) 

After resampling, every sample has an equal length in the direction 
along the scan. 

AY - W (B.2) 

N 

From Figure 3.^, the distance dj may be written, in terms of 
sample widths as 

dj = j (W/N) - 1/2(W/N) (B.3) 

Introducing the term h T ~ Z - 7 * , the distance may also be written 

L j c 

in terms of sample angle y as 

dj » h^ tan ny + h^ tan(U^y - ny) (B.h) 

Equating the right hand sides of (B.3) and (B.U) and solving for 
Uj yields 



(B.5) 
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If Equation (B.l) is substituted into (B.5), the final resulting 


egression is 


n + 1. tan” 1 1 f 2J - 1 h, tan ny + 
Y 2N 


2.1 - 1 h tan(Ny 
2N N 


- tan ny 


The solution of this algorithm requires an iterative approach, 
since the value h^ is not known until the algorithm has been solved 
for U^. As a first approximation. 


Ah-r = hi •* hi 


The algorithm is solved, the updated value h^ is found, and the 

algorithm is reapplied for an updated U value. 

J 

As an example, consider a typical aircraft flight in which 
Z c = 5000 ft. (1.52 km), y ~ 6 mr&d, and symmetric scanning for 
which N - 222, n = N/2 = 111. If it further assumed that - Z^ - 
600 ft. (183 m) , then h-j_ = hjj = 41+00 ft. ( 1 . 3 1 * km). For this example 
the value of U will be computed for J = 50, a value which results 
in a near maximum displacement (see Figure 3.3) and it is further 
assumed that Z^q - 600 ft. (183 m) in the unprocessed array and 
that the ground slope ( X) in the vicinity is 0.5, a relatively 
extreme case. Using these values in Equation (B.6), the first 
solution for Uj results in 



Uj = 111 + 1/.006 tan -1 {(l/l+l+00) [(99 AW (41+00) tan(0.666) 

+ ( 99 /l+l+U ) ( HUOO ) tan( 0 .666) - l+l+OO tan 0.666]} 

Uj = 1+2.6 


For the next iteration the value of Z r = Z^ 2 is utilized instead 




of Z^q. The size of a ground element in the vicinity for this esse 
(from Figure 1.19) is approximately 30 ft. (9*1 m) , and would result 
in a change in elevation over the eight elements of 120 ft. (36.6 m) 
If the ground slope is positive then Z^ 2 = 720 ft. (219.5 m) , 
Utilizing this updated value in Equation (B.6) results in 


Uj = 1+0.9 


a change of approximately 1.7 elements. If one further iteration 
is carried out using Z L = Z^q, the value becomes 

J 


U* = 1+0.7 


If the ground slope is negative, then for the second iteration 
Z^2 = 480 ft. (ll+6.3 m) . Utilizing this updated value in (B.6) 
results in 


Uj = 1+1.8 

and a further iteration produces negligible change. It was found 
during the investigation that a single iteration is sufficient to 
refine the resampled position within one array element. 
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